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Supersonic flow past a bluff body with a detached shock 
Part I Two-dimensional body 


By W. CHESTER 
Department of Mathematics, University of Bristol 


(Received 14 March 1956) 


SUMMARY 

The flow at high Mach number past a body with a rounded nose 
is considered. Viscosity and heat conduction are neglected, 
and the body is assumed to be two-dimensional and symmetrical 
about an axis parallel to the incident stream. 

An exact solution is first obtained in the case y +1, M—> «, 
where y is the adiabatic index and M is the Mach number. ‘This 
solution is then used as the basis of a double expansion in 
5 =(y—1)/(y+1) and M~, after the exact solution has been 
modified to make the series expansion converge near the body. 
The expansion is carried out as far as the terms of order (6+ M-*)?. 

The results are displayed for various values of 6 and M~*; 
typical results are as follows. With M?=0,5=% (y= 1-4), 
and for a parabolic body having unit radius of curvature at the nose, 
the shock is approximately a parabola with radius of curvature 
1-822 at the nose. ‘The distance between the body and the shock 
along the axis of symmetry is 0-323, and the height of the sonic 
point from this axis is 0-744 both on the shock and on the body. 
The actual pressure distribution on the body is shown in figure 4, 
and agrees well with experiment. The pressure falls to zero at 
a distance 0-86 downstream from the nose of the body, measured 
along the axis of symmetry. On the assumption that the pressure 
remains negligible beyond this point, the total drag is 1-39 p,V?, 
where py is the density and V is the velocity of the incident stream. 


INTRODUCTION 

It is well-known that a shock wave is produced when a supersonic stream 
impinges upon a stationary obstacle, and that if the obstacle has a blunt 
nose the shock wave is curved and lies upstream of the body. It is the 
purpose of this paper to investigate the flow behind such a curved shock and 
its position relative to the obstacle producing it. Viscosity and heat conduc- 
tion are neglected (except, of course, at the shock), and the body is assumed 
to be two-dimensional and symmetrical about an axis in the direction of the 
incident stream. 
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An exact solution is first obtained for the case y > 1, M-> o, where y is 
the adiabatic index and M is the Mach number of the incident stream. 
The result for the pressure was originally obtained by Busemann (1933; 
see also Ivey, Klunker & Bowen 1948) using a less direct approach. It has, 
however, been suggested by Stocker (1955) that this is only a rough approxi- 
mation to the truth, even for very large Mach numbers, because the pressure 
field is particularly sensitive to the value of y. 

Since in practice y — 1 is fairly small, particularly at the higher tempera- 
tures associated with transition through a strong shock wave, it would seem 
natural to consider the possibility of improving the exact solution by using 
it as the first term of a series solution in powers of y—1, and in inverse 
powers of the Mach number. _ It is, however, not difficult to show that such 
a solution does not exist, or at least lacks the desirable quality of being 
uniformly convergent in the neighbourhood of the body. ‘The reason is 
that in this region the exact solution fails to give the correct approximation 
(for non-zero 6+ M~?) to the magnitude of the velocity, though, as will 
appear, the direction of the velocity and the pressure are given correctly 
to a first approximation. ‘The first problem is thus to correct the basic 
solution. When this is done, an iterativé procedure can be used which, 
in theory, will give the answer to any desired degree of accuracy. In this 
paper the solution is taken as far as those terms whose magnitude is of the 
second order of small quantities. 

‘The success of the method is due in part to the use of the von Mises 
transformation (Goldstein 1938) in which the stream function is used as 
one of the independent variables. ‘This transformation simplifies the 
equations of motion, particularly when entropy variations must be taken into 
account. It was first used in problems of compressible flow by von Karman 
& Tsien (1938), and has since been used by several writers. It is particularly 
suited to the problem considered here. 

‘The shape of the shock, rather than the shape of the body producing it, 
is assumed to be known a priori. The boundary conditions resulting from 
the transition relations can then be obtained explicitly and applied on a 
known boundary. In principle it is quite possible, and not more difficult, 
to solve the problem for a given body shape; but in practice the former 
approach gives a series which seems to converge more satisfactorily. ‘This 
is probably due to the fact that, if the shape of the shock is known, all the 
boundary conditions, and the boundary on which they are applied, are 
known exactly. When the shape of the body is given, one of the boundary 
conditions must be applied at the body in the region where the error, even in 
the improved approximate solution, is relatively large. ‘This will have an 
effect on the deduced shape of the shock (no longer known exactly), and 
hence on the boundary conditions applied at the shock in the next approxi- 
mation. 

Working with a known shock imposes no serious limitation on the 
results. Most of the analysis is in fact carried through assuming a parabolic 
profile for the shock, and attention is confined to the flow in the important 
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region near the nose. In this region the shape of the body is approximately 
parabolic, the error being within the accuracy of the present solution. 

On the other hand the solution fails to hold far downstream, for it 
predicts that the pressure on the body eventually becomes zero and changes 
sign. ‘This is because one of the assumptions made in obtaining the 
approximate solution is not valid in a region where the pressure, and hence 
the density, falls to a small fraction of its value at the nose. (Essentially, 
the reason is that the expansion of p” in powers of 5 is not uniformly valid 
as p>, and this limits the validity of equation (13).) Experimental 
evidence does suggest, however, that at high Mach numbers the pressure 
eventually becomes negligible. Moreover, the conditions at the nose of 
the body are independent of conditions sufficiently far downstream where 
supersonic flow is well established. More precisely, the subsonic region 
at the nose is independent of the shape of the body downstream of the 
characteristic through the sonic point on the shock. ‘Thus the ultimate 
breakdown of the solution should not affect its validity near the nose. 


NOTATION AND FUNDAMENTAL EQUATIONS 

Let a two-dimensional inviscid uniform stream, with supersonic velocity 
parallel to the x-axis, impinge upon a stationary blunt nosed body sym- 
metrical about the x-axis. ‘The origin of co-ordinates (x, y) is taken at the 
vertex of the shock wave which will be formed upstream of the body. 

The velocity, pressure and density in the uniform stream are denoted 
by (V,9), Po, po respectively. Corresponding quantities in the region of 
disturbed flow behind the shock are (Vu, Vv), pypV*p and pop(y + 1)/(y— 1). 
The variables u, v, p and p so defined are thus non-dimensional. The 
reason for the factor (y+ 1)/(y—1) in the representation of the density will 
appear later. It ensures that p remains finite when y = 1 and the Mach 
number of the uniform stream tends to infinity. 

The equations of conservation of mass, momentum and entropy are 


A _ 


Cc 


} oO . 
a, (pu) + x, (pe) = 9, 


ou ou 0 
Bo Oa = =. 
OX cy pox (1) 
CU oz 60 [ 
ne He ne, 
Ox ry p oy 
© (pe) +0= (pp) = 0 
a= r )-o— i ct 
er ° ay Pp , | 


where 6 has been written for (y—1)/(y +1). 

The first of equations (1) implies the existence of a stream function ¢% 
such that dy, = pu, d, = —pv. In the uniform stream, % = y and is 
continuous across the shock wave. Hence, on the shock, ¢=y. Along 
the x-axis and on the body, % = 0. 


2A2 
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If we change to independent variables (¥, y), equations (1) transform to 


ra Ou 
p _" 


u2+02+(1+8)p/o = 14+(1—8)82M-, 


~——_—_Y 





o /u o/1 
os _ + eee == | 
Oy (;) . ==) m | 

pe’ = f(), 


where f is an arbitrary function of yw. 

It is also suitable to record here the form of equations (1) whenz (= ¢%/y) 
and y are used as independent variables. In the (4, z)-plane the body will 
be z = 0, the shock z = 1, the axis of symmetry y = 0; and we have 


é é ) 
(uy) = —(p+uz), | 
cy Os | 
9 > me re | 
u?+v*+(1+5)p/p = 14+(1-—8)69M-, | 
J 
( P 6 /dz—pu ©) 
oy\pv/ oz\ pu /’ 
! 
pe” = f(y3) J 
Furthermore, 
54 mu — bs 
bee ey am By (4) 
pv pt 
ass ’ ot thet 0 (oy 
Since pv is an odd function of y, it follows that —({—} = 0 on y = 0, and 
; Cy \pv 


hence, by (3), that pu=45z. ‘Together with the second and fourth of 
equations (3) (with f(0) determined by the transition relations across a 
normal shock) this serves to determine all the dependent variables as functions 
of z along the axis of symmetry. 

The distance between the body and the shock along y = 0 is, from (4), 

a ie y 
x =6] lim (2) dz. (5) 

Joy >o0 \p2 
Let the shock be defined by x = X(y) (with X(0) = 0), and let X,(y) 
denote the derivative of X(y). The transition relations across the shock 
wave then give the following boundary conditions to be satisfied at ¢ = y: 


u= 1—(1-—58): — yi ey ||) Gals 


a 
II 
_— 
— 
| 
~~ 
— 
nN 
pd 
A 
— 
ie 
ro 
- 


aad (6) 








Pp = 1 as (1 —$)s7M~(1 +X?) : 
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It follows paeinem that 
1 5 M-2 ) ‘a \ (1+ 6)/(1-—4) 


fd) = 1-3)4 Teyep ~ Tos | --(1-8)3-4-M-*(1+X3 oC) f 
| (7) 





From the third equations of (2), we now have 
1 


eles 
II | 
~ ~ 
| a =— 
Pea 2| 
Oy. SO, 
B= Ble 
——— ead 
+ + 
rr, ~~, 
Cis <I 
| 
Blo 
“aa” 
° 


With the help of equations ae this gives 

u 0 (#1 
- = X,(y)+d5— | —d 
= = X,0)+85, [sab 


“yp 


(8) 


N 
candies, 


: . ~&a dz 
= X,(y)+8 | = dz+—. 
ro Jz0Y pv pv 
It follows that the equation of the eet is 
x = X(y)+8 € — dy. (9) 
J g pv 
Note that, as y + 0, this becomes 


Pe ie y 
) | lim (=) dz, 
~Oy->0 pv 


which by (5), is just the distance between the body and the shock. 


‘THE BASIC SOLUTION 
We consider first the case 5-> 0, M-?- 0 (such that 6-1M~ remains 
finite). Equation (7) then becomes 
1 
Ih; = —___ +6 1\f-2: 
f(!) Tt XA) ; 


and, when this is combined with equations (2), we arrive at the following 


relations: 
é é Sei Xp 0 [u 
Sot. Hien Se, J <8 (10) 
oy = os 14+ X2(f)’ ob\v 
Hence u/v = X,(¥), which also follows immediately from (8). ‘The shock 
wave thus coincides with the body in the physical plane (though not in the 
(4, y)-plane). There is, however, a discontinuity in the rest of the flow 
variables between the shock and the body, in addition to the discontinuity 
across the shock. In the (,y)-plane, the solutions of the first two of 
equations (10) give 














7 X,() “ y — X,() 
{1+ XG(p) 91 + XH(y) pr?” {1+ XE) PPL + XG(y) 7?” | (11) 
1 df X, i Xb)dp 
P= TFX3Q) dy (4X9 Vj T+ X2(p)}!2" 
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In particular, on the body (4% = 0), we have 
Toe fe 
p = sin?8B+sinB— | cosBdy 
dy J0 


where § is the inclination of the tangent to the horizontal. ‘This is the 
result originally obtained by Busemann (1933). 

Henceforth the shock is assumed to be.a parabola with unit radius of 
curvature at the vertex, so that its equation is x = }y*, and equations (11) 


give 
py Be. ob 


(1 + p2)2(] +. y2)t2? (1 + f2)2(1 + y2)12” 
(1+y2)12 (1+ 42)32 | 


J ote 
P< (eye? PO yA TP) J 








u= 


(12) 





In order to improve this approximation, the straightforward approach would 
be to substitute this basic solution in the neglected terms of (2) and solve 
the new set of equations as far as the terms of order 6 and M-*. Unfor- 
tunately this process does not converge near /=(. For example, as 
ys» 0 v ~ & according to (11); and hence the third of equations (2) would 
give u/v ~ dlogys. ‘The reason is that the solution is not expansible as a 
power series in 6 and M-? when { = 0. It will appear that the velocity 
components are in fact O(8 + M-*)!? on the body. ‘Thus the basic solution 
(12), which gives u = v = 0, is not a valid approximation to u and wv for 
ys = 0, in the sense that the neglected terms are small compared with the 
term retained. On the other hand, p, p and u/(vy) are all finite and non-zero 
on the body according to (12); and hence for these quantities the basic 
solution is presumably a valid approximation in the above sense, for the 
error is presumably o(1) when (6+ M-?) is small. This will in fact be 
verified when the higher order terms are obtained. 

Equations (12) are accordingly used as the basic solution for p, p and 
u/(vy). It is then a simple matter to improve the corresponding solutions 
for the velocity components. 

From equations (2), (7) and (12) we have, correct to the first order, 


E = flv) p'-* = f(d){1 +28 log p} 


] ‘ 
= ——— +§61M-2(1 + J?) —6 —-2M-2(14+- pf?) + 
ernie 51 M—~(1 )-—o 1~*( ib?) 
ee | ee 
+5,;0+ M (1 + ?)} loe( 75) | ’ (13) 
so that 

u2+y? = ja J? + 3{5+ M-*(1 +7)! log A (14) 

. 1+ ia dns 7 "\ 1442 P 


When this is combined with u = yz, and all but the leading terms are 
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discarded, the following uniformly valid approximations are obtained : 








y ome nl? 
u= (1+y?)2(1 + pp |” + 3d log(1+ 5 ) 
: [ ) 1/2 t (15) 
a < f+ 3d 1+y2) , | 
(1 +y?)! 2(1 +?) 2 \¥ ¢ log( +) ) | 


where d has been written for (6 + M-*). 


FIRST-ORDER APPROXIMATION 
From equations (12) and (15), we now have 
5 (d+ M-* *)(1+y?)* ; (16) 
pu = {? + 3d log(1 +7) 9(1 + ¥?)’ 
and this is uniformly correct as far as terms which are O(d). ‘This is 
substituted in the right-hand side of the equation 


dfu\ a8 
ob\v}/  —s dy \pv’ 


which is then integrated to obtain the next approximation for u/v. The 
boundary condition is, from (6) with X,, = 4, 


5 M-* 
2 =yt 5 (I +?) + —(1 +")? (17) 





9 


C 
The details are straightforward and are omitted; the only care required is in 
ensuring that the terms which are discarded are all o(d). 
The result is 


a dipy(1 
vo? {f? + 3dlog(1+y 


+ y?) _ (14 2\o (18) 
yWElog(1 ty) 20 tI)» 





where 


ub 4,2 
= . a i 2 3dloo(1 42) 4 
g = 4dsinh 3dlog(1 +52)" 2d log 3dlog(1 +?) 
g 5 g a 


d y? 1+¥y? 
[ee sae, a iy > —2 y - 77 2 C 
‘2 raat +) 1] 2(d—M )log| m Z| i (i 


the term in u/v which is singular at y = 0. We have already seen that, near 
y = 0, pu ~ dz = dys/y, or 


u= U+ 








diy(1 +?) 2 
(1 + ys?) ?log(1 + y?) 
say, where U is O(y?). Equations (15), (18) and (20) then combine to give 


U 
—=y —y(1+y")g; (21) 





(20) 


and the singular term which appeared in equation (18) has now disappeared. 
These relations can now be substituted in equation (14) to give 
Bua = diy re 
TPAC PYF log +95) * 


a5 x1 =f) er, ie 1+y? 1/2 
(1 SRY sical poe A) | ; (22) 
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and hence, from (20), 
dipy 
Ty og +95) 
y(1 -g) 1+ y?\"]!2 
. M~ 23 
+ aH X1 + giz yi? + 3(d+ M- ¥*) logl te » (23) 
— diy? 
(1 + y?)"?(1 + f?)!? log(1+y?) ~ 
(1+gy?) ; mi? 1+y2\12 
f z eae ae -2, /,2 y Z 2 
(1 ry)! a1 +P) 3 yb 3(d T M Ww )log(; + ya ° (24) 
‘The quantities of prime physical interest are now easily obtained. The 
shape of the body is calculated from (9) using (16). The equation is referred 
to an origin at the vertex of the shock, and so automatically gives the distance 
between the body and the shock. ‘The equation is, correct to the first order, 
4y? ) 
3d(1 + y®) log (1+y%){ 
+43M-*(1+y?)?log(1+y?). (25) 
Finally, the pressure is calculated from the equation 
op ou 


Cus — oy . 


u= 








o 
c 











x = Ly2+41d(1+ y*)? log: 


the right-hand side being given by (23). When this is integrated and 
combined with the boundary condition p = (1 —54)/(1+y?) for % = y, one 
finds that 

(1+ y2)32p (1+ 2)l2 


; : us 
— d{? log(1+ y2) + 3v?4+4—- t af2)1/2] sinh-} = 
= d{5log(1 + y?) + 3y? + 4—4(1 + f?)!?] sinh E dog aa] 


j 2 : jj2 2 
- df} log(1 +4) +4y*+ 2} log 5 ey at td ]+ 











2 3d log(1 +?) 1+ (1 +y?)!?| 
i eee ee ae 
+. + yf2)h2) 2— 7 r 
ah+%) | 2-9 il: 2log tsa 1+y*) log(1+ ah} 
+ d(1 +y2)!%(y-?—3) —}[ fl? + 3d log( 1 +32)}12— Ye 
jloglt+y)  ¢ 
at = 
2d | aren ainhé 7 


1+ - ¢ 9\3/9 
+ M-*(1 + f?)"2 | 847 7y? + jlog(7 7, =) | —7M~*(1+47)8?. (26) 


SECOND-ORDER APPROXIMATION 


In principle the approximation can now be carried as far as the second 
order terms by repeating the process with the first-order approximation 
replacing the basic solution. 

Some care is required in the calculation of the integrand that occurs in 
equation (9) for the body contour. An expression for v is required which is 
correct as far as terms which are O(d3") on the body. Equation (24) is not 
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accurate to this order, since terms have been omitted from within the square 
brackets which are O(d?), but which do not tend to zero as {-> 0. On the 
body, their contribution to v is O(d*?). Essentially, this means that the 
first-order solution must be corrected in the same way that the zero-order 
solution is corrected before proceeding to higher order terms. ‘Thus in the 
second of equations (2) we use the first-order solution for p to get an improved 
expression for u?+v*, analogous to equation (14). This is then combined 
with the first order solution for u/v to obtain improved expressions for the 
velocity components. 

In practice the algebra becomes prohibitive once the stage of calculating 
the pressure is reached. Since this investigation is mainly concerned 
with conditions in the neighbourhood of the nose of the body (which are 
independent of the behaviour far downstream), the calculation is accordingly 
restricted to the expansion of the second-order terms omitting fourth and 
higher powers of y; for this purpose equations (3) have been found superior 
to equations (2). Actually, the complete expression for the second-order 
terms is obtained in the case of the equation of the body contour, and is 
quoted in the appendix ; in equation (27) below all fourth and higher powers 
of y are neglected. Calculations suggest there is little loss of accuracy in 
ignoring the higher powers of y at least as far as the sonic point on the body. 
The final results only are quoted below. 

The equation of the body is 


4 1 4 3 4 4 
sae eS hime. ws bw tie, 4 ea 
x = ddlog zd +d ( log 3d 3) }dM~ log aq +3M—4 (icg id 1)+ 


4 ., 5, nf17, 54,7, 4 161 
al iE + dlog —_ ld+ IM +a°( Fog 3a * 16 °8 39 ~ 6 == 


3 
-ant(Sog#, -2)+01-(ioe—2)]+009. @n 


In the following expression for the pressure on the body ( = 0), the first- 
order terms are quoted exactly; in the second-order terms fourth and 


higher powers of y are omitted. Thus, 


ee 3 phlog(i+y"*) é : 
(1+y?)?*#p-1= al 5 | sane" ((1 +y?)!2—1}{y-?-3}-1- 
a 
— 62 2 3 —u 
{2 log(1 +?) +597} log ¢ sleet 5} ; log(1 +7) + 


+ 2{? log(1 +?) + $y?+ +2} log b+ (1+ 32)" | mn 
+ M-*[8 + 7y?+3log(1 +y%)—7(1 +9?)°?] -3M-*°(d— M~) + 3a? — 


(21, 4 «15. 4 43) 
oe: 2 2 2 
%y | 4 | 16!°8 3a * 76 !°8 3a + 96/7 


(3 & 7) .3 f 4.) 
M-*. — log =- M+; log = —4; |. 
+dM 5 0854 + mi *4 1 08 37 i} ] (28) 
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NUMERICAL CALCULATIONS 








The results embodied in figures 1 to 4 are all based on equations (27) and 
(28). It was found that the equation of the body was given with sufficient 
accuracy by equation (27). For example, the neglect of the higher-order 

T T = T oa 
| | 
| | | | 
| | 

| | Moi XLay | 
| | M =O-17d x | 
| . -2 
\ on ee | M’=027d 
| ee a "| 
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Figure 1. The radius of curvature of the shock on the axis of symmetry for various 


values of d (=(y 
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The unit of length is the radius of 
curvature of the body on the axis of symmetry. 


The isolated points are 





experimental values. 


terms in y? produces an error in the ordinate of the sonic point of only 1°% for 
y = 14, M*=0. On the other hand, the pressure was calculated from 
equation (28) as it stands. 
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Figure 2. The distance between the body and the shock along the axis of symmetry 
for various values of d (=(y—1)/(y +1)~M-*). The unit of length is the 
radius of curvature of the body on the axis of symmetry. The isolated points 
are experimental values. 
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The previous analysis is based on the assumption of a parabolic shock with 
unit radius of curvature at the vertex. In practice the results are of more 
significance when displayed for a given body, and in figures 1 to 4 the results 
have been adjusted so as to apply to a body which is a parabolic cylinder with 
unit radius of curvature at the vertex. 

The radius of curvature of the shock, when y = 0 is given in figure 1. If 
this is denoted by 7, the equation of the shock would then be y? = 2rx. 
Since the adiabatic index decreases at the high temperatures associated with 
the conditions behind a strong shock wave, it was thought to be of interest to 
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Figure 3. The ordinate of the sonic point —————— on the body -—---- on the 
shock. When M-? = d, the sonic point on the shock is at the vertex. 


exhibit the results for various values of 8. On the other hand, since the more 
significant parameter is d (= 5+ M~-*), and since 0 < M~ < d, the variation 
with dis shown only for the two extreme values of M-?. 

Figure 2 shows the distance between the body and the shock, along the 
line of symmetry, as a fraction of the body radius of curvature. Again the 
variation is plotted against d for the two extreme values of M~. 

The ordinate of the sonic point is shown in figure 3 for both the body and 
the shock. On the body the pressure at the sonic point is given by 

p, = (1+8)-(1 —8)(1+ 8-6 M-?)-8- 9/8 (29) 
This result can be deduced from the second of equations (2). Equations (28) 
and (29) are together used to obtain figure 3. ‘The ordinate of the sonic 
point on the shock follows at once from (6) with XY, = y/r. Explicitly it is 
§(1 — M-2)] 12, . 
r| 1+5M | 
When M-? = d(or 8 = 0), the sonic point is at the vertex of the shock. It is 
clear from this fact and from figure 3 that the position of the sonic point on 
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the shock (for a given body shape) is far more sensitive to variations in both 
dand M~ than the position on the body. When M-? = 0, for example, the 
variation on the body is less than 20°, of its minimum value over the range 
considered in figure 3. In this connection it is interesting to study figure 
4, which shows the pressure distribution along the body when M-? = 0, 
and y = lor 1-4. ‘The former is the exact solution (1+ y?)-3?; the latter is 
calculated from (28). Figure 4 shows that the pressure in the neighbourhood 
of the sonic point is even less sensitive to variations in d than the position of 
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Figure +. The pressure distribution on the body; the abscissa is the distance along 
the axis of symmetry, measured from the stagnation point, as a fraction of 
the body radius of curvatures ————— M-?=0, d= 1/6 (y = 1-4). 
-- M~* = 0, d= 0 (y = 1). The crosses denote experimental values 
for a circular cylinder with M = 5:8. 

the sonic point itself. Indeed, the Busemann solution would seem to be 
rather better than has previously been suggested. Moreover the maximum 
error between the stagnation point and the sonic point (and for some distance 
beyond) occurs at the stagnation point, where the pressure is known exactly 

from (7) and the second of equations (2). Explicitly it is given by 
Ps = (14+8)7 {1-8 —8(1 —8) M2} 40-9), (30) 

The experimental results displayed in the figures are taken from a paper 
by Oliver (1956). In addition, the results of an experiment by Holder 
(private communication) are shown in figures 1 and 2. Both experiments 
were carried out on circular cylinders at Mach numbers of 5-8 (Oliver) and 
4 (Holder). A value of 1-4 for the adiabatic index has been assumed in 
plotting these results. 

The theory tends to underestimate the stand-off distance and the radius 
of the shock, judging by the available experimental evidence. The pressure 
distribution on the body near the nose, however, agrees well with the 
measurements made by Oliver on a circular cylinder. 
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When d = M-* = 0, the pressure on the body is always positive and 

approaches zero asymptotically at large distances downstream. For non- 

zero values of the parameters, however, the pressure becomes zero at a 

finite distance downstream; and this is borne out by experiment. On the 

assumption that beyond this point the pressure remains negligible, the total 
drag on the body would then be (including the upper and lower halves) 


2poh *| p dy, 


the integral being taken from the stagnation point to the point of zero pressure. 

When d = M~ = 0 the total drag is just 2p)V’? for a parabola with unit 
radius of curvature at the nose. The corresponding result for d = 6, 
M~=0 is 1:39p,V?. The pressure becomes zero when y = 1-31, or 
x—X,) = 0-86 (here x—x, is the distance downstream measured along the 
axis of symmetry from an origin at the nose of the body). 


APPENDIX 
The equation of the body contour, with the second order terms given 
exactly, is as follows: 
x = hy?+ $dP(QO-—R)+3M?P!R+ ia PtQ*sR + 7y*?) +3M4*PYO-1- 


~(1—y2)R?—(7 + 6y°)R + 14P4P—1)} + 4d? P 4 in 1—4y°)R+6P— 
—2-—2(P+1)?-—(3R+ by? + 8)S—2y?R+-3 j * Ecosech dé | + 


+4d?P# | + 2y?)R? —8+4{1—3y?—-y?R1-—2(P—1)R4}(P+1)?- 
—3(4y? —-5)R+ 4y(1 — 3”) PS + (S?4+2S—R+ 2R1S)(3R + 6y? + 8) + 
“aR 
+ | {6+2R +3 log(4ty* coth 3€)}é cosech € at | - 
~O 


—1dMP*| | (7 — 4y2)R + 16 + 12y?— 14 P?}0 — (2y-2 + 24+ 18y2)R— 
—4(1—y2)R? + 2(3R + 6y? +8+14P)S+ 4(1+4y?+7y!)(P+1)7- 
—2R-1( P?— P34+4R)— 6f “expe —1)geoseché a], 


4y2 1 
where P=(1+y*)", Q= log, ae 


R = log(1+y2), S = log }{1+(1+2)#2}. 
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On the theory of hypersonic flow past plane and 
axially symmetric bluff bodies 
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SUMMARY 
The ‘ Newtonian-plus-centrifugal’ approximate solution 
(Busemann (1933) and Ivey (1948)) for hypersonic flow past 
plane and axially symmetric bluff bodies in gases with the ratio 
of the specific heats y constant and equal to unity is rederived 
using ‘boundary layer’ techniques together with the von Mises 
variables x and %. A method of successive approximations then 


‘ 


gives a closer approximation to this solution for e = (y—1)/(y +1) 
small and the free-stream Mach number infinite. Formulae for 
the streamlines, shock shape and pressure distribution are deter- 
mined to this approximation. ‘These formulae are valid for any 
plane or axially symmetric shape, giving the ‘stand-off’ distance 
of the shock wave from the body as }elog(4/3e) and « times the 
nose radius of curvature for plane and axially-symmetric flows 
respectively. Particular results are computed for a number of 
special shapes. For certain shapes, the theory has a singular 
point where the first approximation to the pressure vanishes 
(9 = 60° for a sphere). Actually, the theory is not applicable 
where the pressure becomes too small. The corresponding 
theory for gases of general thermodynamic properties is deduced, 
the approximation being valid provided the total energy of the 
gas is large compared with the energy contained in the trans- 
lational modes of the gas molecules. 


1. INTRODUCTION 

The term ‘hypersonic’ has been coined to describe flow regimes where 
the free-stream Mach number is considerably in excess of unity. Theoreti- 
cally, several useful approximate solutions have been obtained for uniform 
hypersonic flow over slender bodies where the bow shock wave is attached 
and the fluid is assumed to be inviscid (Van Dyke 1954, Lees 1956). For 
bluff bodies, however, the bow shock wave is detached, and a mixed 
subsonic-supersonic region exists behind it. ‘This makes the theoretical 
solution of the problem considerably more difficult. Nevertheless, by 
making the additional assumption that the ratio of the specific heats y of 
the gas is near one, an approximation to the pressure on the surface of a 
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symmetrical bluff body in an inviscid flow can be obtained (Ivey, Klunker & 
Bowen 1948, Busemann 1933). ‘This is known as the ‘ Newtonian-plus- 
centrifugal ’ solution for reasons which will be seen below. 

Behind an oblique shock wave in a uniform flow of velocity U, density p, 
and Mach number M, the normal component of velocity w,, is given (§ 2) by 


u, = r— U sin + O( 55) 





y M? 


/ 
as [M -> oo, where the gas is assumed to be perfect with constant specific 
heats, and where © is the angle the shock makes with the oncoming flow. 
We see that as y > 1, uw,» 0. Thus, in the limiting case of a strong shock 
wave and y near unity, it is possible for the bow shock wave to rest actually 
on the body, since the boundary condition for an inviscid fluid, that the 
normal component of velocity is zero, is satisfied. Fluid will still cross the 
shock, however; for the density behind the shock is 


P = Pod a a (sn) , 
which becomes infinite as y > 1, and hence the mass flow puw,, across the shock 
wave remains finite. This fluid can be squeezed between the shock and the 
body, since the stream-tube area which is proportional to pyU/pg becomes 
small for y ~ 1 (q is the velocity behind the shock wave and remains finite 
as y-> 1 due to continuity of the tangential velocity u,= U cos across 
the shock wave). ‘The pressure behind the shock wave is 


» 
p= Po U® sin’? ~ p,U? sin?0 
y+ 


as y->1. This was given as the pressure on the body surface by Epstein 
(1931); but, as we shall see, it is incorrect for bodies of finite curvature. 
The name ‘ Newtonian ’ is given to this result, since it is the pressure on a 
body surface placed in a stream of inelastic particles as given by Newton 
(1689), the particles losing all their momentum normal to the surface on 
impact. It would seem that to attempt to imply a closer relationship is 
questionable, for the highly compressed gas discussed here bears little 
resemblance to the “‘medium rarum quod ex particularis quam minimis 
guiescentibus aequilibus et ad aequalis ab invicem distantes libere dispositis 
constat”’ which Newton postulates. 

For bodies of finite curvature, however, the pressure on the body differs 
from that upon the shock, since there is a finite amount of fluid in the narrow 
region between the two. The pressure difference is equal to the centrifugal 
force on this layer of fluid. ‘Taking this into account leads to the ‘ Newtonian- 
plus-centrifugal’ solution given by Ivey, Klunker & Bowen (1948) and 
anticipated by Busemann (1933). 

In this paper, we will follow these authors in using the ‘strong shock’ 
approximation for large M, but will assume ¢« = (y—1)/(y+1) to be small 
but not zero, with the solution for « = 0 as a first approximation in the 
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solution for infinite Mach number. ‘This successive approximation 
procedure is something like expanding in powers of €; the solution will 
be obtained as far as the first term in e. 

In §2, the conditions across the shock wave are derived from the 
Rankine—Hugoniot relations. A further simplification is introduced by 
assuming the gas to be perfect with constant specific heats. In practice, 
of course, the conditions for such an assumption to be valid will not prevail. 
However, the theory will be developed throughout using this assumption 
in order to obtain a definite result. ‘The modifications to be introduced 
if this is not so are discussed in § 7, where the solution for a gas with arbitrary 
thermodynamic properties is given. 

The ‘ Newtonian-plus-centrifugal’ solution for two- and _ three- 
dimensional flow is then obtained in §3 from the equations of momentum, 
continuity and energy for a compressible inviscid flow without heat 
conduction. We assume that between the shock wave and the body there 
is a thin ‘boundary layer’ in which changes perpendicular to the body 
surface are large compared with those along the body. As for other 
boundary layers, we assume the layer thickness is small compared with 
the body dimensions and that the velocity component normal to the body 
is small compared with that along the body, or, alternatively, that the 
streamlines deviate only slightly from the body shape. A measure of the 
magnitude of the variables being given by their values on the shock wave, 
the density will be of the order of py/e and the pressure of the order of p,U. 
We must not, however, be led into making too close a parallel between 
this boundary layer and the classical viscous boundary layer where in most 
cases any centrifugal effects are negligible. Introducing the ordinary 
boundary layer coordinates (x, y), the momentum equation in the x-direction, 
approximated in the above manner, states that the velocity component u 
in the x-direction is constant along streamlines. (This is because with 
e small the enthalpy changes very little in an adiabatic expansion.) ‘The 
momentum equation in the y-direction gives the pressure gradient normal 
-to the body as due to the centrifugal forces produced by body curvature 
alone. From this, the pressure at any point of the layer is obtained as a 
function of the stream function ¢ (or, in axially symmetric flow, the Stokes 
stream function) and the coordinate x. Using the energy equation in the 
form of constancy of entropy along streamlines, the density can also be 
obtained in terms of these variables. ‘lo this approximation the shock is 
assumed to have the body shape. Finally, the position of each streamline 
can be derived by integration from a knowledge of p and wu as functions 
of x and ys. By replacing y by its first approximation on the shock wave 
in this expression, we obtain the shock shape. ‘This result is given in §4 
and plotted in figure 3 for the cases of the circular cylinder, parabolic 
cylinder, sphere and paraboloid of revolution. Knowledge of the stream- 
lines of the flow and the shock shape enables us to calculate the small 
component of velocity v parallel to y and a higher approximation to u, 
and hence obtain a second approximation to the pressure ($6). In 
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particular, the variation of the pressure on the surface of a sphere is 
computed and plotted in figure 5. 

The velocity tangential to the shock wave increases as we move away 
from the stagnation point; and hence, on many bluff bodies (those on 
which the curvature does not decrease with distance from the stagnation 
point), the centrifugal force on the fluid will increase with increasing 
distance from the stagnation point. ‘The pressure rise across the shock wave, 
however, being proportional to the normal velocity, will decrease. Con- 
sequently, the pressure drop between the shock and the body may eventually 
annul the pressure rise across the shock wave, and the pressure as calculated 
from the theory will fall to zero on the body. ‘The point where this happens 
is a singular point of the theory, the approximation not being valid there. 
This difficulty arises from the failure of the first approximation at this 
point. It is no longer possible to assume that the stream-tube area is still 
small enough for the shock wave to rest approximately on the body surface 
once the pressure on the body has fallen to a small fraction of pyU?. The 
problem then arises of finding the first approximation to the shock shape 
in this region. It would seem impossible to do this without making some 
further assumption, since, immediately we admit the layer between the 
shock wave and the body to be no longer thin, we are beset with a further 
difficulty of having two velocity components of the same order of magnitude. 
The form which this further assumption should take remains a matter for 
conjecture. ‘The streamlines as they emerge from the narrow layer 
presumably fan out to cover the region between the shock and the body. 
Several assumptions have been tried as to the nature of this fanning-out, 
such as allowing the gas to separate from the surface, but none seems to 
give results of the magnitude required from experimental evidence. ‘They 
will not, therefore, be further discussed. ‘The restrictions implied by this 
singularity can best be seen by consideration of a few examples. Un- 
fortunately, it would appear that the circular cylinder and sphere are 
particularly bad from this point of view, the approximation breaking down 
at 54°44’ and 60° respectively from the front stagnation point. On a 
parabolic cylinder and paraboloid, the position is much more favourable, 
the singular point being situated at infinity in both cases. However, even 
on these shapes a position is reached where the pressure falls to a value 
too small for the theory to be valid. 

The assumptions underlying the theory introduced in this paper, 
namely, that the Mach number is infinite and the ratio of the specific heats 
of the gas is nearly unity, are only approximately true in practice. In the 
hypersonic regime the ratio of the specific heats for air reaches a minimum 
of about 1-2. At large Mach numbers (and hence high temperatures) 
dissociation of the gas is almost certain to be well advanced. Variation of 
the molecular weight of the gas with temperature and pressure and, more 
important, the large variation in specific heats must then be taken into 
account. This theory can be extended (§ 7) to cover gases of quite arbitrary 
thermodynamic properties, provided the parameter ¢ is replaced by 
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K-! = py/p, (where p, is the density on the shock wave). ‘The parameter K 
is not necessarily constant, however. ‘The solution obtained in §7 holds 
for large values of K. 

An interesting feature of the flow discussed here is the apparently large 
shear in the layer between the shock and the body, even on the assumption 
of inviscid flow. ‘This shear is almost linear (figure 4). Moreover, the 
component of velocity u, since it is constant to the first approximation 
along streamlines, is zero upon the body itself (as the body surface is the 
streamline through the stagnation point). ‘his being so, the non-slip 
condition at the surface for a viscous fluid is already satisfied to this 
approximation. It would seem that, on this account, viscosity may not 
play such an important part in the flow as expected. ‘The theory does, 
however, neglect heat conduction, which will be especially important in 
regions near the stagnation point where gas temperatures are large. 
Nevertheless, this may not have too disastrous an effect, as we may expect 
heat conduction to lead to lower temperatures at the body surface and 
hence a thinning of the boundary layer. 

A further difficulty which arises in considering a real gas is that the 
equilibrium assumed behind the shock wave may take some time to achieve. 
Relaxation times may be quite large for vibration and dissociation, which 
would be excited at the high temperatures the gas is likely to experience. 
This would cause the ratio of the specific heats to vary considerably in the 
flow immediately behind the shock wave. ‘The relaxation time being 
roughly proportional to the density, the importance of this effect is 
determined by the size of the parameter (pd) !, where d is a length 
determining the scale of the flow. 

The results given here seem to agree fairly well with a solution given by 
Lighthill (1956) in cases where A 1s large. ‘This solution refers to conditions 
near the stagnation point of a sphere, comparison being made only in that 
region. ‘The ‘stand-off’ distance « of the shock wave on a sphere disagrees, 
however, with the value Se given by Schwartz & Eckermann (1955). ‘This 
may be due to their making some assumptions with regard to the flow near 
the stagnation point. ‘Their experimental results would seem to confirm 
neither value. Experimental results given by Oliver (1956) (on flow past 
a cylinder with a hemispherical nose) for the pressure on the body surface 
seem to agree more closely with the original ‘ Newtonian’ solution than the 
modified one. ‘The Mach number in these experiments was, however, 
still quite low (about six). 

During the preparation of this paper for publication, the author’s 
attention was drawn to the (then unpublished) work by Chester (1956). 
Chester also treats the present problem by a process of successive 
approximation, starting with the ‘ Newtonian-plus-centrifugal’ solution. 
It was decided, however, to publish the two results separately as his approach 
to the problem is different from the present one, and also the further 
development of the theory proceeds along different lines. In so far as it 
is possible to compare Chester’s results with those given here, they are 
identical, 
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2. ‘THE SHOCK CONDITIONS 
An oblique shock wave makes an angle 0 with a uniform flow of velocity 
U’, the density, pressure and enthalpy in the flow being pp, py and % 
respectively (figure 1). Behind the shock the velocities tangential and 
normal to the shock wave are u, and u,; and the pressure, density and 
enthalpy are p,, p, and 7, respectively. Normal to the shock wave the 
Rankine—Hugoniot equations (stating continuity of mass flow, momentum 


Figure 1. The variables at the shock wave. 


and energy) hold. ‘Tangential to the shock wave the component of velocity 
is continuous. From these relations, together with a knowledge of the 
thermodynamics of the gas, we obtain the conditions behind the shock 
wave. ‘Thus, 
ppUsin® = p,u,, 
Pot pol? sin? = p, +p, Re, (2 1) 
io t+ SU? sin?O = 1,4 }u' 
Ucos® = u,, 
and ty = 1( Pr» Ps), (2.2) 
the latter equation being a purely thermodynamical relationship. For 
hypersonic flow the dynamic terms on the left-hand side are much larger 
than the state variables p, and i. Thus, the equations (2.1) may be 
approximated as 
ppUsin® = p,4,, 
pul? sin? = py +p, U;, (23 
roy 2 2 3) 
4U?sin?O = 1,4 4u%, 


UcosO = u,. 


Ww 


B2 
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Alternatively, equation (2.3) can be written in terms of K = p,/py as 


u 


Usi 


and 
Also, 


In general, the shock equations (2.4) are solved for K by substituting /, 
and p, into the thermodynamic relation in the form p = Q(p,7). If K is 
large, however, a method of successive approximation is found to converge 


very 


for given 7, the ratio 7, p,/p, (which by (2.5) determines A) varies only 
slowly with p, or with p, and so can be determined by trying successively 
different values of p, and p,. This ratio is 3/2 times the ratio of the total 
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u, = Ucos0. 
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rapidly; for 7, is known from (2.4) to be very nearly }U*sin?@, and 





Ipy to the translational energy of the molecules, which would be 


expected to be large for real gases in hypersonic flow. Hence, the 
4 S 


assun 


iption that A is large can be interpreted as ip/p being large in the 


region behind the shock wave. 


For K large, we see from equation (2.4) that the normal component of 
velocity is small, and the density large behind the shock wave. 


For a perfect gas with constant specific heats, AK is a constant independent 


ot (-), 


and equals « ! = (y+1)/(y—1). ‘Thus the assumption of large K 


/ 


requires y to be nearly unity for such a gas. ‘The theory will be developed 


for this ideal gas in the following sections. 


3. "THE EQUATIONS OF MOTION 


Consider flow past a symmetrical bluff body either plane or axially 
symmetric, placed in a uniform stream of velocity U. ‘The shock wave 


S 


F) B 


ay Ab 


Figure 2, The curvilinear coordinate system used in the boundary layer. 


S shock; B, body. | 
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in front of the body will be detached and lie a distance 6—usually referred 
to as the ‘stand-off’ distance—upstream of the stagnation point. We 
introduce the usual boundary layer coordinates (x, y) (figure 2). The y-axis 
is normal to the body surface, and the x-axis lies along the body in the 
plane formed by the normal and the direction of the uniform stream. 
A further coordinate s is introduced on the body surface, such that (x, y, =) 
form a right-handed system of axes. ‘The elements of length in the (x, y, 3) 
directions are denoted by hdx, dy and kdz respectively. The case of two- 
dimensional flow can be reproduced by putting k = 1. 

The equations of momentum, continuity and energy may then be 











written u ou Ou 1 a 
- = += +uvK+ — — = 0, 
h ex ev hp ex 
U2 o@ lo 
= eo Pa ee Fh. = 0, | 
h ox eh p oy [ 
' z (3.1) 
2 
ao (Apu) +r — (hkpv) = Q, 
cov 
uoas es 
7 a- +05— = J, 
h ox ey 
where S is the entropy and « the curvature of the lines y = constant. We 


assume now that the stream tube area is small behind the shock wave when 
e is small, and it is therefore possible, in the limit « +0, for the shock to 
rest on the body. For small « we assume a priori that the distance of the 
shock from the body is O(e), and that slope of the shock relative to the body 
is O(e). (In the two-dimensional case it will be seen that it is later necessary 
to modify this assumption. ‘These magnitudes then become O(¢ log ).) 
On the shock, u,/U = O(e) and u, U = O(1); and thus, resolving in the 
direction of the coordinate axes, we have u/U = O(1) and v/U = O(e) on 
the shock, also p/p,U2 = O(1) and p/p) = O(c) from (2.4). In the flow 
behind the shock wave, these variables will have the same orders of 
magnitude. Changes across the thin boundary layer between the shock 
and the body will be large compared with those along it. Assuming that 
the above conditions hold, except perhaps in regions near the stagnation 


point, we introduce the new variables 


f p : €p F u , Uv 
p = —> p=-, uu= —, v= —, (3.2) 
pol x Po l Ue 
and P v va se 
: d’ _ ed ' 


and then the dashed variables are all O(1). Substituting in (3.1) and noting 
that A,k and « are O(1), we obtain 





u’ ou’ ou’ 

— — +v’— = Ofe), 

h cx’ oy’ (<) 
l o 7 
— Pu = O(€) (3.3) 
p’ oy 
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Assuming the gas to be perfect, we have S = p/p’ and hence 


u’ do ([p’ , 3 (p’ ; 
—_—_— — _ + WU om — a ( y " 4. 
h se(5) sA5) Ie log €) (3.4) 


Neglecting terms of higher order than the first, we therefore have, in terms 
of the original variables, 





uct ou 

. vx = U, 

h Cx oy 
1 < ee 
- == UK (3.5) 
poy 


a Ox\ p ov\ p : 


with the continuity equation remaining unchanged. ‘The first equation 
of (3.5) states that the velocity component uw is constant along streamlines, 
and the final equation states that the enthalpy is constant along streamlines. 
‘This is consistent with Bernoulli’s equation which is, to this approximation, 


w+(5") = 1U?, (3.6) 


The second equation of (3.5) demands that the pressure gradient normal 
to the body shall be that due to the centrifugal force produced by the 
curvature of the body alone. 

Introducing a stream function » to satisfy the continuity equation, we 
require that 


ous Cr , " 
kpu = —, hkpv = — fl (3.7) 
ay a 
Equations (3.5) then state that u and p/p are functions of ys alone, and that 
Cp uK ; 
2 nl 8 
ae (3.8) 


4. THE EQUATION OF THE SHOCK AND STREAMLINES IN TWO AND THREE 
DIMENSIONS 
The mass of fluid flowing across the shock wave is pu, where u,, is the 
velocity normal to the shock. Now, 


pu,, = p(ul+vm), (4.1) 


where /, m are the direction cosines of the normal to the shock wave. By 
using (3.8), this may be written 
; ; 
1/135 ous 
pu, = + ~— +m,—}, (4.2) 

k\h ox ey 
where /,, m, are the direction cosines of the tangent to the shock wave. 
If s is the distance measured along the shock wave, therefore, 
1 dbs 
k ds* 
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On the shock wave pu, = p»Usin©, where © is the angle of the shock 
to the oncoming flow, or, in terms of s, 


dn 
= pU—, 4.4 
pu,, Po a ( ) 
where 7,(x) is the distance of the shock from the line of symmetry. Hence, 
from (4.3), dy, 1 dibs - 
a = (4.5) 


ds k ds’ 
or us = poU | kdy.. (4.6) 


Let us now assume % = 0 on the line of symmetry; and then, since k = 1 
in two dimensions and k = », in three dimensions, we have 
wu = p,Un,, 
” ts a (4.7) 
2P0U%)s > 
in two and three dimensions respectively. If (x) is the distance of a point 
on the body itself away from the line of symmetry, then 


», = (x) + Y(x)cos M(x), (4.8) 


where © is the angle of the body to the uniform flow direction and y = U(x) 
is the position of the shock. Thus, if ‘/(x) is known, equation (4.7) together 
with equation (4.5) give % as a function of x on the shock wave. On the 
assumption that !/(x) is small, however, 
ub = pol (x) 
= bpUnXx) (4.9) 
in two and three dimensions respectively, if we neglect O(e). Thus we 
know & explicitly on the shock wave: ys = pyU V(x), say. 
Also, © = P+tan{h-\(d//dx)}, 
which gives on neglecting O(e), 
© = ®. (4.10) 
(hus the shock conditions (§2) reduce to 
u = U[cos®+ O(e)], ) 


v = O(Ue), 
.... *. | (4.11) 
Po € 


P_ = sin + Ole), 


pol : 
which together with (4.9) give the dependent variables as functions of # 
on the shock wave. 
If we now introduce a new variable € defined by % = p,UY(€), then we 
can replace ys by € which is the x coordinate of the point where the stream- 
line ys crosses the shock wave. 
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From the equations of motion, u and p/p are functions of %, and hence 
of €, alone. ‘Thus, using (4.11) we obtain 


u = Ucos®(€), (4.12) 
and ea eU*sin®@(€) 

p 
throughout the flow field. 

Also, from (3.8) 
p=|7 4 
; (4.13) 
"uk 

or p = pAx)+ | = dys, 


where the suffix s refers to conditions on the shock wave at a station w. 
Now, « =: «,, 9+ O(e) and k =k, _,+O(e); and hence 
p = p,+ Z | u dis 
gaol (4.14) 
p + pol | ut’’(t) dt, 
a 
where u = u(x, = pyUV(t)) and «, 4 =k, etc. From (4.11) and (4.12) 
we obtain 
Pp. = po? sin’D(x) and u= Ucos®(€); 
and therefore 
P = = sin*@(x) + ms ‘t''(t)cos M(t) dt. (4.15) 
pol : k, ‘ 
This is the result given by Busemann (1933) and Ivey (1948). Again 
using (4.12), we obtain 


ZL 


| sin2D(x) + |” ‘Y'’(t)cos P(t) dt 
pP ty - 


Peta sin2@(é) 


It is now possible to reintroduce the coordinate y into the problem, 


(4.16) 


z 





For, from (3.7), we have 


dis, (4.17) 


y= 7 
- | » Rou 


keeping x constant; or, making the previous approximations and introducing 
(4.16), we obtain 
Y'(t)sin2@(t) dt 





*/ sin2@(x) 4 2 7 ‘t'’(s)cos D(s) dS + cos P(t) 


2 


for the equation of the streamlines. 

The above integral does not converge for the two-dimensional case, 
however, since cos ®(t) = O(t)as t--0. In the three-dimensional case the 
integrand has no such singularity, since the function ‘Y’ = O(t) also. 
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In the two-dimensional problem, therefore, it is necessary to obtain a 
better approximation to u which will be valid for ¢ or € small. It will be 
seen that magnitudes of the order of eloge will replace those taken to be 
O(e) in the previous approximation; and thus the error terms previously 
denoted by O(e) must be replaced by O(eloge). Bernoulli’s equation in 
the form (3.6) may still be expected to hold near the body, since v = 0 by 
the boundary condition there. ‘The energy equation (3.1) for a perfect gas 
with constant specific heats may be written 


== (4.18) 


where ps and px are the values of p and p where the streamline crosses the 
shock wave. Equations (4.18) and (3.6) together then give 


‘ l+e px Pp ry ‘ g 
ju + = | ag 1+ 2elog-—} = $U*[1 + O(elog*e)]. (4.19) 
: 2 /} Po \ Px, y 
Also, the pressure on the shock wave is, from (2.4), 
we 1 dy) 
p = (1—e)pol alba + tant _ , (4.20) 


Substituting in (4.19), we obtain 


uy? = U| cost0(e cali loe( Oe 





) +24'(E)cot O(E) > + 


- O(elog*s) | , (4.21) 


where py = p(€,&) and p is given by (4.15). 

The term in curly brackets is O(eloge), and becomes important near 
the body surface where the first term is small. Since, however, this higher 
approximation is important only near the body surface, it will be sufficient 
to substitute its value actually on the body to obtain the position of the 
streamlines from (4.17). Thus, putting 4 = € = 0 in the second term 
of (4.21), the required approximation for the velocity component u is 


° 1/2 
a= U | cost) — 2e log | : (4.22) 
py U? 





where p(x, €) is given by (4.15) (since /'(€) ~ Oand p(€, €) ~ py U? as € -- 0). 
In the two-dimensional case this will replace the original value of u in (4.17). 
Thus the equation for the streamlines € = constant is given by 
sin°@(t) dt 





ii | » R(x, t){cos?@(t) — 2e log R(x, 0)}!2’ daa 
mnie R(x, t) = sin*(x) +’ {” cos ®(s)sin @(s) dS 
in two dimensions, and ‘ 
r> n(t)sin®@(t) dt (4.24) 


ala |, O(x, t)cos O(t)’ 


mane O(x, t) = sin?@(x) + = ("cos (s)sin D(s)n(s) ds 
n Je 
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in three dimensions, the variables ky, x, and 7’ having been replaced by 
n(x), —®' and sin® respectively. On the shock wave, € = x, and hence 
the equation of the shock wave is obtained by replacing € by x in (4.23) 
or (4.24). 

The validity of three formulae near the stagnation point is doubtful 
since the approximation that one velocity component is much larger than 
the other may not be true there, for both components of velocity are small 
in the neighbourhood. It will be shown later ($5), however, that the 
formulae remain valid. Using (4.23) and (4.24), therefore, we can obtain 
the ‘ stand-off’ distances in two and three dimensions as 

6 | 


+ 3 
= ele. (4.25) 
3 


a € 


and a € (4.26) 


by noting that »(¢) ~~ ¢ and ® ~ }7—14 a for ¢ small, where a is the radius 
of curvature of the nose of the body. 

Geometrically, the simplest types of blutf bodies are, of course, the 
circular cylinder and the sphere. For these, ® = 57—4 and 7(x) = a sin8, 
where 6(= x/a) is the angle in the plane (or meridian plane) measured 
from the front stagnation and a is the radius. ‘The equations (4.23) and 
(4.24) then give 


r—a ! cost dt 
2¢ — re 
a » [3cos?4— cos*t] {sin — 2 log( 1 —$sin?@)}!* 
é f/ —- 2sin?(&/a) aoe sin*(€/a) 
ee 7 a % log eT a eee 3 cos*4 log | r 3 as a 
(3 cos? — 1) ellog(1 — ? sin)! 2—3sin°é 
(4.27) 
to this approximation for a circular cylinder, and 
r-a ' cos*t dt 
€ | - ——; 
a Jy 4(sin 36+ sin*t) 
é 4+ 8 (24? 4 5 Be 
: | (1 —*)log —z-— Hlog{ | _ 
7 v1 He ~ OF 


, By 3 
\ Stan | (4.28) 
%—f 
) 


for a sphere, where x = ¥ (sin 34), 8 = sin(€/a) and r is a polar coordinate 
measured from the centre of the sphere or cylinder. The shock shape is 
deduced by putting € a = 4 in (4.27) or (4.28). 

‘These formulae become useless near the singular point of the theory, 
given by R(aé,0) = 0 for (4.23) and O(a, 0) = 0 for (4.24). For a circular 
cylinder this is ¢ = 54 44’, and for a sphere 6= 60. For a parabolic 
cylinder or paraboloid of revolution, however, this singular point is removed 
to infinity. Equations (4.23) and (4.24) then become 


\ | y- - ie dt 
2h c( = yy (+P e+) '4 3 log(1 + ¥? 46%)|"? 
\ 


(E246) log[l + (¥2/482)]/ siti 
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for a parabola, and 
’ 2\ 3/2 a/b 

Re = 2¢( + a3) “a en. Eee (4.30) 

2b 4b?) J, (14+2?)!?[G.(¥/2b) + G_(d)| 
for a paraboloid of revolution, where G__(t) = ¢(1+#?)'? + sinh”! t, (X, ¥) 
is a system of Cartesian coordinates (with XY along the line of symmetry 
and Y perpendicular to it) in which the parabola has an equation Y* = 4bX, 
and = is the value of Y at the point where the streamline ¢ crosses the shock 


Thus the shock wave is obtained by putting = = Y in (4.29) or 


wave. 
(4.30). 
These results are plotted in figure 3. 
af 











| 1 
oO | | | ke, 
os 5 tO {-> 2-0 
Figure 3. The distance “/ of the shock wave away from the body as a function of 


the distance x from the stagnation point for (a) a circular cylinder (4 (£1) 5 
(b) a parabolic cylinder (y — 1:1), (¢) a sphere, and (d) a paraboloid of 


revolution. a is the radius of curvature of the body at the stagnation point. 
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Figure 4+. The first approximation to the velocity profiles in the boundary layer on 
a sphere at @ = 20° and 30°. 
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5. BEHAVIOUR OF SOLUTION NEAR THE STAGNATION POINT 

As has been pointed out in the previous section, the formulae given 
there may not be valid near the stagnation point where both components 
of velocity are small. ‘That these formulae do in fact remain true can be 
shown in the following manner. Near the stagnation point it is possible 
to consider the flow as incompressible, since the velocity components will 
all be small there. The simplification in the previous theory occurs because 
we are able to approximate the velocity term in Bernoulli’s equation by 
including only one component of velocity. It is now no longer possible to 
do this. Nevertheless it will be seen that a somewhat fortuitous piece of 
cancelling enables us to overcome this difficulty. 

Consider first the two-dimensional problem. ‘The momentum equation 
in the normal direction is 


u, Cu lop u; 


~ 
| 
| 
| 
— 
ut 
ea 
— 


where (7,6) are polar coordinates with origin at the centre of curvature of 
the body near the stagnation point, and p is assumed constant. In the 
previous work we neglected the two ‘convective’ terms in the equation. 
Now, however, the velocity component wu, is no longer small, and thus we 
approximate the equation as 


or, introducing the stream function ¢/, 
ee C2 
a7 (P+ 2pu;)— — 0, (9.3) 


where i is chosen so that 
Oo c 
— tax 
cr p Je uw 


It will be noticed that the pressure in equation (3.8) has now been replaced 
by p+ spuz. Consequently, if we consider Bernoulli's equation for in- 
compressible flow, viz. 

pt ip(u;+usz) = pq +P,, (5.4) 
where py and gy are values on the shock wave, we see that replacing the 
pressure by p+ Spu> reduces the problem once more to one containing the 
single component of velocity u,. ‘Thus, integrating equation (5.3), we have 


. i — 

pt+tpur—| — db = p.+ spur (5.5) 

. vy, 7 ; 
: & 

where the suffix s denotes values on the shock wave at a station 6. 

Substituting for p+ pu? from equation (5.4) into equation (5.5) and 


differentiating with respect to 4%, we obtain 


9 Cly 
iors 
ee) ous 


O es Uy a 
ay (2% re =: (5.6) 


U 
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Hence, by the cancellation of the term }pu?, the problem is reduced once 
more to one containing u, only. Approximating the shock wave in this 
neighbourhood by X¥ = —c+(¥?/2d)+O(Y*), where (X, Y) are Cartesian 
coordinates with origin at the stagnation point and the X-axis along the 
line of symmetry, we then have 

gz, = &U2+(U2Y2/d?) + O(Y’), 
and 


Pe = (1—€)p,U2[1 —(¥2/d2)] + O(Y?). (5.7) 


On the shock wave, 4 ~ p,UY; and hence, retaining only first order in & 


in (5.6), we obtain 


u, a = ab — Bu,, (5.8) 
where 
% Jn -2e+0/(e)] 
pid? 
2 Feith 
i o()), 


and a is the radius of curvature of the body at the stagnation point. Solving 
(5.8), we then obtain 


p* Biles — Mes (5.9) 


’ 


ue x —Bu,—u? 
neglecting terms which are o(1) on the right-hand side. Approximating 
(5.9) further, we have 
uu— —— = 3€U?sin?0, (5.10) 
assuming that 
d = a{1+o(I)]. 


Also, since 


(r—a)/a = | | [pu,]-} db, 


eh lelog( s—aeraraa) (5.11) 
a 


3€p7 U?a? sin? 


then 


which is the limiting form of (4.23). 

If we proceed in a similar manner for the three-dimensional case, we 
confirm the result (4.24) in the region of the stagnation point, but with 
an error term which is O(e*?) and not O(¢*log’e) as result (4.24) indicates. 


6. 'THE SECOND APPROXIMATION TO THE PRESSURE 
In §4, we obtained a second approximation to the shock shape and 
streamlines. From this result we can hence obtain the values of the 
variables in the region behind the shock wave to the second approximation. 
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To do this, we must consider a higher approximation to the equations of 
motion (3.5). “The momentum equation in the y-direction taken from (3.1) 
and written in terms of the variables (x, us) becomes 


u C2 C 
— —u*x +uR 8 0), 
h ox Cus 
or Cp uk 1 oz 
I (6.1) 


Oy k hk ox 
In §3, the second term on the right-hand side was neglected completely. 
We must now seek a higher approximation to the first term and a first 
approximation to the second term. Integrating (6.1), we obtain 


p = p.(x)+ | ‘ (+ — 2 =) dis. (6.2) 
From §4, we know the equation of the streamlines, 
y = f(x, §), or F(x, i), (6.3) 
say. ‘Then the equation of the shock is given by 
y = ee) = Wiz). (6.4) 


‘Thus the angle the shock makes to the oncoming flow is {® +tan71(Y(x)/A)}, 
whence the pressure on the shock is 
Pp = py U1 —€)sin?/® + tan-(/'(x)/A)} 





: se Yf'(x) 
p, Us (1 — €)sin*® + sin 20 —— (6.5) 
h,(x) 
to the next approximation. 
From Bernoulli’s equation, we have 
P (yu? + v?) U2, 
y l p - = 
— ve, 
ueniana u> U? -( vo tale); (6.6) 
€ p 
also, from the energy equation, 
Pp _ pe 
p Px 
where px, px denote values at the shock wave. ‘Thus, 
u? U2 a3 § te 2elog® ) +0(¢) 
Po px 
ae x, xX ; Y'(x. 
U*< cos?@(v,) — 2e sin?D(x log A ) — sin 20(x,) | hie’ 7 
{ . Px(x,) hy(x.) 
(6.7) 


where x = x, is the coordinate of the point on the shock where the streamline 


crosses it, and is related to the stream function by 
ys = py U[n(x,) + f(x, )cosP)x,)] 
= JpyUln(x,) + (x,)oos (x)? (6.8) 








of 
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in two and three dimensions respectively. Equation (6.7) then gives the 
second approximation to the velocity u. 
From equation (6.3), it follows that 





cus C F C F 
a eae yy? 
CX OX] Ow 
Or 
1 oF joF 


(6.9) 


phk ox | Ob’ 

Since we only require a first approximation to 7, we can replace phk in (6.9) 

by its first approximation, and hence obtain 

e Py, oF /oF 

hi, Ra Do Pex Ou 

where P denotes the first approximation to the pressure given by (4.15), 
and P, its value on the shock wave where the streamline %& crosses it. 

Let us now introduce a variable € to replace the stream function %, 


+ o (Ue), (6.10) 


putting 
is = py 7(€) 
» py Un?(E) (6.11) 
in two and three dimensions respectively, and let us write this ys = pyU‘Y(€) 
as before. ‘The equation (6.10) becomes 
eU P«/cf /ef\d¥ : 
—~ (2 2 — +0(Ue). (6.12) 
h,k, P \ox/ 0& 

‘The variable x, is then related to € by 


n(E) = (x) + U(x, )cos (x,), 


N(E)cos M(E) 


\ 6 = ee, (6.13) 
© 
and equation (6.7) can be w ritten 
: W'(E) 
u U cos @(€)< 1 + Y(E)@'(E) — > tan VE) — 
h(€) 


2 


—€ tan*(2)loe( 5) to(e)>. (6.14) 


Noting in equation (6.2) that us, corresponds to € = x +[/(x)cos ®(x)/'(x)], 
we then obtain the pressure to the second approximation as 








({’ K t/ uu 
i sin’? — e sin? + he sin 2) — ke 7 = cos?® + 
"ST wx, ) ; var 2 (C)tan D(C) 
1 ‘os D 1+ Y(0)0'(c =o eet 2 
Ea (¢) (C)®'(¢) hid) 


© Mig hy 6x\ Pall) hy Ro ex] 66) de 





{P C C 2 ea /of of i) 
— etan2(Z)log P| ‘ P(x,¢) 1 (af /ef\d" i 


(P55 i 
(6.15) 
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where «(x,€) = «(x,y = f(x, 0)) = Ko(x)[1 — Ky f(x, 0] to this approximation, 
k= 1 or 7(x)+f(x, O)cos ® and 


Pp sin? — (@’/kp) | ds) ds 
Py rT? sin*p(C) 

The arguments of the functions in (6.15) have been left out when they 
are simply x. ‘The error in (6.15) is O(e?log*e) for two-dimensional flow, 
and O(e??) for three-dimensional. 

The expression (6.15) can be evaluated for any body shape when the 
function @(xv) is known, Inthe case of asphere, ® = }7— 6, and the pressure 
on the body surface becomes 





p 








—~. = Pit Po, (6.16) 
pol i 
where sin 34 
Pi = 35ind’ 
and 
Ps = cos?6{ 2A’ tan 6-1} + A(6)cos 20 — 2(1 — } sin?@) - 
1 1 / sin 36+ sin°¢ 
— | A(d)cos3¢ddh+ - | cos*é log] ———- = |d¢+ 
sind !,, sind /, 3 sin 4 cos*d 
9 . ay (/ sin*y(sin®d — sin®ys) dys 
ee + SIT) Of — ly ae 
sin 4 | (sin 34+ sin®ys)? 
oag | Cos*ih(sin3é — sin3ys) dys 
+ 2.cos?364 | —— ee FET ae 
» (sin 34+ sin®,)3 
with "  cos*d dd 
\(4) | 


' sin 38+ sin? © 
The final three integralsin p, and the function A can each be obtained in terms 
of simple functions, but lack of space prevents us from quoting them in 
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Figure 5. The pressure coefficients p, and p, on the body surface, and the shock wave 
for a sphere. ‘Total pressure p = poU*(p,+ep,). S, shock; B, body. 
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full. The functions p, and p, are plotted in figure 5. We see that near the 
singular point where p, is zero the pressure coefficient p, has a non-integrable 
singularity, i.e. pp ~ (sin 3@)-83 as @-+ 47. This formulae can therefore not 
be expected to hold above 6 ~ 35°. It is also possible, if a little care is 
exercised in dealing with the stream function, to obtain from (6.16) a simple 
expression for the pressure along the axis of symmetry, which is 


cee 4 


giving a check on the result. 





7. GASES WITH ARBITRARY THERMODYNAMICAL PROPERTIES 


Under the conditions in which the previous theory holds, it becomes 
impossible to treat the gases as perfect with constant specific heats. It is 
necessary therefore to generalize the preceding theory to include gases 
which have quite arbitrary thermodynamical properties. In fact, it will 
be found possible to deduce the previous results for such a gas provided 
that we have a knowledge of how the enthalpy varies with pressure and 
density. 

The shock equations can be written in the form given in equations (2.4). 
Provided, therefore, we know 7 = i(p,p), these equations can be solved to 
give u,, u,, p and p on the shock wave. For K large, the magnitude of 
these variables is obtained from (2.4), and we can therefore assume as 
before that the shock wave rests on the body to a first approximation. 
Hence we obtain u and wv on the shock wave. Behind the shock wave the 
variables will have magnitudes of this order and, in particular, p/zp is small. 
We make the same approximations as before in the ‘ boundary layer’ between 
shock and body. ‘the final equation of motion (3.1), the energy equation, 
written in terms of the enthalpy and in the coordinates x and %, is 
a anil (7.1) 

OX pox 
the derivatives being taken along streamlines. Since p/ip is small behind 
the shock wave, equation (7.1) becomes, to the first approximation, 


© 


t 
_* 0, (7.2) 
or 1 = I(€), say. 

Again, the velocity component normal to the body is small, and Bernoulli’s 
equation states 


i+4u2 = 102, (7.3) 


Thus wu is constant along streamlines. ‘The momentum equation in the 


y-direction remains unchanged, and hence the pressure is given by (4.15) 


as before. Thus we know p = p(x,&) and i = I(€) from the shock 
conditions and (4.5). The enthalpy of the gas being given as a function 
of pressure - density, we thus have the density behind the shock as 

= Q[p(x, €), 1(€)], say. Approximating the stream function as before 


F.M. 2C 
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in (4.9), and using the integral (4.17), we obtain the equation of the stream- 
lines. In the two-dimensional case, however, it is again necessary to obtain 
a better approximation for u near the body, i.e. for € small. This is obtained 
from (7.1) by not completely neglecting the second term on the left-hand 
side, but replacing it by its first approximation. 

‘Then, 


i=1@+( 2a, (7.4) 





where p(x, ) = Q[p(x, £), 1(Q)], and p is given by (4.15). Equation (7.3) 
then gives the required approximation to uw. 
The formulae replacing (4.23) and (4.24) can then be written 
% sin?(t) dt 
“OU Xx, t = J OR 1/2’ 
| cos*@(t)—2 | oi on 
Po t Jqp(%,s) ox J 
where p(x, s) = Q(p)U?R(x, s), I(s)), with R(x,s) as defined in (4.23), for 
the two-dimensional case, and 





(7.5) 








Ca > 7(t)sin®@(t) dt 


| 9 p(x, t)cos D(t) 4 (7.6) 


where p(x, t) = Q(p)U?Q(x, t), 1(t)), with Q(x,t) as defined in (4.24) for 
the three-dimensional case. The formulae (7.6) and (7.7) then give the 
equations of the streamlines € = constant. On putting & = x, they give 
the equation of the shock wave. 

Having obtained the shock shape y = ‘/(x), it is then not very difficult 
to generalize the expression for the pressure (6.15) in a like manner to 
obtain 


Ky(x) U(x) d¥ 


—=— >— cos*®(x) + 


= sin?®(x) + H(M(x)) + ca sin 2(x) — 











py U? h(x) ko(x) n(x) dx 
on yr 
+| li cos (0), 1+ (00) — 2 tan o(e)— 


sec’*@(Z) (7 1 OP 
U2}. p(x, t) dx 








eae ee es aad 
hy Ry 0x | p(x, C)Mo Ro \ax/ ag dt tla dl, (7.7) 


where y = f(x,€) is the equation of the streamlines from (7.5) or (7.6), 


p(x, t) = Q(P(x, t), Z(t)) with P defined in (4.15), and 


K(x, ¢) = K(x,y = f(x, f)), 
etc. asin (6.15). The pressure on a shock wave of angle © has been written 
Pp = p,U*[sin?O + H(O)] (7.8) 


to this approximation, being deduced from the shock equations (2.4). 








— «-=— A wet te oP OP 








), 





On the theory of hypersonic flow past bluff bodies 387 


8. CONCLUSION 

The theory developed above would seem to be useful in predicting the 
flow over many bluff bodies, both plane and axially symmetric. On all 
bodies, it gives a solution in regions near the nose for flows in which the 
Mach number is large and the density ratio across the shock wave large. 
However, difficulties arise in regions where the pressure falls below a certain 
value of the order of p,U?/K. In these regions the solution is no longer 
valid. ‘To obtain a solution at points beyond such a region, it would be 
necessary to make further assumptions about the flow. 


The author is greatly indebted to Professor M. J. Lighthill for much 
valuable help and encouragement throughout the preparation of this paper. 
The author is also grateful to the Department of Scientific and Industrial 
Research for a grant during the period of this research. 
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SUMMARY 


‘This note advocates a model of the steady flow abouta bluff body 
at large Reynolds number which is different from the classical 
free-streamline model of Helmholtz and Kirchhoff. It is suggested 
that, although the free-streamline model may be a proper solution 
of the Navier-Stokes equation with » = 0, it is unlikely to be the 
limit, as 4 —> 0, of the solution describing the steady flow due to 
the presence of a bluff body in an otherwise uniform stream. The 
limit solution proposed here is one which gives a closed wake. 

A closed wake contains a standing eddy, or eddies, whose 
general features can be inferred from the results of an earlier 
investigation of steady flow in a closed region at large Reynolds 
number. In all cases, the drag (coefficient) on the body tends 
to zero as the Reynolds number tends to infinity. ‘The procedure 
for finding the details of the closed wake behind two-dimensional 
and axisymmetrical bodies is described, although no particular 
case has yet been worked out. 


1. INTRODUCTION 


‘The determination of the steady flow about a bluff body placed in a 
uniform stream of incompressible viscous fluid at large Reynolds number 
(where the word ‘steady’ in this context implies that somehow turbulence 
has been suppressed) is an old problem, for which no completely satisfactory 
solution is available. Enough is known of solutions of the Navier-Stokes 
equation for us to anticipate that at large Reynolds number viscous forces 
will be small (when made non-dimensional in the usual way) everywhere 
except in the neighbourhood of a number of singular surfaces of the fluid. 
Away from these singular surfaces, or shear layers, the flow is essentially 
inviscid, and various analytical and numerical techniques are applicable; 
near the singular surfaces the flow has a boundary layer character, and again 
there are many relevant aids to analysis. However, these principles and 
techniques are of no avail unless the properties of the singular surfaces (viz. 
their position, and the velocity increments across them) in the fluid are 
known, and this is the principal stumbling-block of the problem. 

In the case of fluid streaming steadily past a body of such streamline 
form that no separation of the boundary layer on the body surface occurs, 
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the properties of the singular surfaces are evident enough, and a complete 
self-consistent description of the flow can be given. ‘The body surface is 
one singular surface, and another is the surface formed from all the stream- 
lines that extend downstream from the trailing edge in the corresponding 
irrotational flow. In the limit of infinite Reynolds number, the velocity of 
the fluid is everywhere the same as it would be for an inviscid fluid, except 
at the body surface and on the wake surface, and the difference between the 
solution for » -> 0 and that for » = 0 (where p is the viscosity of the fluid) is 
here trivial. 

In the case of bluff bodies immersed in a uniform stream, and in a number 
of other kinds of flow in which boundary-layer separation occurs, the 
properties of the singular surfaces are far from being evident intuitively, and 
cannot be calculated directly even for the simplest body shapes. Nor does 
experiment provide much guidance, since the steady flow in the wake of a 
bluff body is unstable, and becomes turbulent, at a Reynolds number which 
is not large enough to reveal clearly the character of the asymptotic steady 
flow. It is necessary, therefore, to make some guess about the properties 
of the singular surfaces before the character of the flow at large Reynolds 
number can be determined. ‘The problem has no direct physical impor- 
tance, of course, in view of the ever-present instability, but it has some 
mathematical interest for those concerned with properties of solutions of 
the Navier-Stokes equation. Moreover, it may well be that a knowledge 
of the steady flow in the limit of infinite Reynolds number would allow the 
determination, by some kind of asymptotic expansion, of the flow at the 
upper end of the range of Reynolds numbers at which the flow is stable, 
more readily than by an expansion valid in the neighbourhood of zero 
Reynolds number. 

The model of the flow that is commonly used is the free-streamline 
model, based on Helmholtz’s (1868) concept of vortex sheets, and worked 
out fully for the case of a (two-dimensional) plate set broadside-on to the 
stream by Kirchhoff (1869) and later by Rayleigh (1876). It is unlikely 
that these authors were concerned to find the complete form of the limit 
flow, and their aim may have been simply to use the observable phenomenon 
of vortex sheets shed from a body in order to provide a more realistic 
representation of the velocity distribution near the plate than can be obtained 
from a wholly irrotational flow. However, be that as it may, the free- 
streamline model of the limit flow is the only one available at the present 
time and is commonly regarded as a correct representation of the flow that 
would occur at large Reynolds number in the absence of turbulence. It 
is in this capacity alone that it is here being criticized and replaced by 
another model. 


2. THE HELMHOLTZ—KIRCHHOFF FREE-STREAMLINE MODEL 
The free-streamline model is undoubtedly right in its division of flow in 
the limit of infinite Reynolds number into regions of inviscid motion 
separated by singular surfaces across which the velocity and its derivatives 
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may be discontinuous. In this respect we can improve on Helmholtz’s and 
Kirchhoff’s notions only by recognizing that vortex sheets need not always 
be uniform. Apart from this feature of vortex sheets, or detached boundary 
layers, which nowadays we should take for granted, the essence of the free- 
streamline model of flow past a bluff body isthat it supposes the velocity of the 
fluid to be zero everywhere inside what may be called a wake bubble, the 
boundary of the wake bubble being a singular surface. ‘Two other pieces 
of infermation must be supplied before the model determines the flow field 
completely; one is the location of the curve (which reduces to two 
points in cases of two-dimensional flow) of intersection of the body surface 
and the surface of the wake bubble, and the other is the (uniform) pressure 
in the wake bubble. It is not clear whether these two pieces of information 
are independent or not. So far as can be judged from work on two- 
dimensional flows, choice of the wake pressure seems to determine the 
position of intersection of the body surface and the wake bubble and vice 
versa; for instance, Southwell & Vaisey (1946) were able to find a definite 
shape of the wake bubble behind a circular cylinder, by a numerical solution 
of the equation for inviscid flow, for each of a number of different assumed 
wake pressures and with no assumption about where the free-streamlines 
met the surface of the body. 

Kirchhoff found that for the two-dimensional flat plate that, when the 
wake pressure is equal to that far upstream, py, say, and the free-streamlines 
spring from the edges of the plate, the wake bubble is open, with a width 
which increases as x'* at a large distance x from the body. Levi—Civita— 
see Birkhoff (1950)—found this same asymptotic parabolic form of the 
wake bubble for a family of two-dimensional bodies with the wake pressure 
taken as py. It appears that, for some choices of the pressure in the wake, 
the bubble is closed at the downstream end. Closed wake bubbles, which 
are necessarily characterized by a cusp at the downstream end in order to 
allow the velocity just outside the bubble boundary to be uniform, have 
been obtained for a circular cylinder (Southwell & Vaisey 1946) and for a 
truncated aerofoil (Lighthill 1949) by choosing wake pressures larger than pp. 

The present objection to the free-streamline model of steady flows in the 
limit of infinite Reynolds number takes one of two different forms according 
to whether the wake bubble is open or closed. If it is open, the objection 
is to precisely that property of the model—that is to say, the model appears 
(to this writer) as wrong in conception. _ If it is closed, the objection is that 
viscous stresses at the boundary of the bubble would build up an appre- 
ciable circulation in the wake bubble and that the model is not self- 
consistent. 

The objection in the latter form rests on the results of a recent investi- 
gation (Batchelor 1956) into the properties of steady flow in a closed region 
at large Reynolds number R. It was argued there that viscous stresses 
acting in the shear layer at the boundary of the closed region have a small 
but persistent effect on the flow in the closed region, and that, provided the 
flow is truly steady, however large the Reynolds number may be (i.e. 
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provided the limit operations t > oo, R-—> o, t being the time measured 
from an instant at which arbitrary initial conditions are prescribed, are 
carried out in that order), the flow in the closed region has a unique form. 
The form of the vorticity distribution could be determined explicitly in the 
cases of closed flows that are two-dimensional or axially symmetrical, 
although a multiplicative constant remained to be determined in each case 
from the condition of matching of the ‘inviscid’ flow in the core of the 
cavity with the viscous flow in the neighbourhood of the boundary. In 
one typical, but mathematically very simple case, this multiplicative 
constant was calculated, and it was established that, as expected, the general 
level of the velocity in the closed region was of the same order of magnitude 
as that of the velocity of the surrounding boundary. ‘The implication of 
this work for truly steady flow past a bluff body at very large Reynolds 
number is that, if the wake bubble is closed, it will contain a standing eddy, 
or eddies, in which the velocity distribution has a unique and (in sufficiently 
simple cases) calculable form, and in which the velocity is of the same order 
of magnitude as (although no doubt smaller than, by a factor of, say, two or 
three) that outside the wake. This motion in the wake bubble will have an 
appreciable effect on the shape of the boundary of the wake bubble, and the 
notion of an equi-pressure boundary to a closed wake is not self-consistent, 
except perhaps as a rough approximation in suitable cases. ‘The properties 
of a closed wake, with allowance for the internal circulation, are considered 
in the next section, this being just the model that is advocated here. 

We turn now to a consideration of the other of the two alternatives, that 
is, the free-streamline model with a wake bubble which is not closed at the 
downstream end. The objections to this picture of the flow about a bluff 
body in the limit of infinite Reynolds number spring mainly from one’s 
distaste for such a drastic interference with conditions at infinity far down- 
stream. I cannot see any reason why a free-streamline flow pattern with 
an open wake should not be a valid solution of the Navier-Stokes equation 
in the limit of infinite Reynolds number (unlike a free-streamline flow with 
a closed wake, which is not self-consistent), but I think there is room for 
doubt about whether it is the solution that corresponds to the boundary 
conditions appropriate to a bluff body placed in an otherwise uniform stream. 

One’s suspicions about this are aroused by the sudden changes that 
appear to be necessary in the general character of the flow when the Reynolds 
number is changed continuously. At sufficiently small Reynolds numbers, 
the streamlines in contact with the two sides of a two-dimensional body 
certainly meet again further downstream, enclosing a finite region of fluid 
behind the body (see, for example, the calculations and observations of the 
flow past a circular cylinder, at Reynolds numbers between 20 and 30, 
described in Goldstein 1938, §20). How then could the region bounded 
by the two branches of the surface streamline change from being closed to 
being open, as the Reynolds number increases? It is not sufficient to 
suppose that the point of closure moves further downstream, and ultimately 
to infinity, as the Reynolds number increases (nor is there de‘inite evidence 
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from calculations that this does happen, once the Reynolds number is above 
some small value, such as about 30 for a circular cylinder), because a closed 
wake bubble, however long, would develop an appreciable internal motion 
by viscous action and the bubble boundary would inevitably take up the 
one possible shape—which (so it is argued in the next section) is such that 
the length of the wake bubble is finite when the Reynolds number is infinite. 

‘he boundary conditions that must be specified in order to make a 
solution of the Navier-Stokes equation unique are not known with certainty, 
so that it is not possible to give a precise statement of the mathematical 
problem under discussion. ‘There is a fair presumption, however, that we 
are looking for the limit, as R —- 0, of that steady solution of the Navier— 
Stokes equation for which the velocity of the fluid is zero at the surface of 
the given body and is uniform everywhere at infinity. It seems to be true 
of all the existing numerical solutions of the Navier-Stokes equation at 
finite Reynolds number that this set of boundary conditions is necessary 
and sufhcient for uniqueness. If this presumption is correct, the free- 
streamline model with an open wake bubble gives an unwanted limit 
solution, in view of its property of non-uniformity of the velocity at infinity. 

Pursuing another line of thought, when the Reynolds number is very 
large, but not infinite, the singular surface at the boundary of the open wake 
bubble is converted into a shear layer of small but finite thickness. ‘The 
thickness of the shear layer will ultimately increase as x!? (in a case of two- 
dimensional flow) by viscous entrainment, like any boundary layer with 
uniform external conditions, and the fluid in the wake bubble must supply 
part of the inflow to the shear layer. A small back-flow (from the region far 
downstream, and towards the body) in the wake bubble is thus necessary, 
and, since the width of the open wake bubble also varies as x"? (in at least 
some cases), this backflow can be thought of as a uniform flow inside the wake 
bubble. ‘The interference with the flow conditions at infinity that is 
implied in the free-streamline model now appears as even more drastic when 
the Reynolds number is not infinite, since the presence of the body must be 
supposed to be responsible for the generation of a small back-flow origi- 
nating far downstream. Is it not likely that the solution represented by 
the free-streamline model with an open wake bubble is that which is 
obtained on/y by the imposition of suitable boundary conditions at infinity, 
i.e. by the external imposition of a small uniform backward velocity (which 
is of order R-!? and hence becomes zero in the limit of infinite Reynolds 
number) over a region y? < const. x x atw—> oo? 


3. THE PROPOSED LIMIT FLOW WITH A CLOSED WAKE 
The discussion of the Helmholtz—Kirchhoff free-streamline model 
given in the preceding section has already indicated the general features of 
the preferred alternative model of limit flow past a bluff body. The point 
of view taken here is that an open wake bubble is inadmissible for the 
problem at hand, although it may be a feature of the limit solution for 
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some different set of boundary conditions at infinity. We should there- 
fore seek a limit solution for which the wake bubble is closed, this being the 
essence of the present proposal. Opinions may vary about the admissi- 
bility of an open wake bubble, but I think there would be general agreement 
that, if a solution describing steady flow about a bluff body in the limit of 
infinite Reynolds number and exhibiting a closed wake bubble could be 
found, that solution would be regarded as the desired solution in preference 
to the free-streamline solution. I have not been able to find such a solution 
mathematically, but I believe that its existence is suggested by a number of 
arguments. 

For the moment, let us accept the suggestion that the desired limit flow 
past a blutf body possesses a closed wake bubble and explore its implications. 
There is an immediate implication that the drag coefficient for the body is 
zero (in the limit, R —» o), since the flow everywhere outside the combined 
body and wake bubble is inviscid and irrotational. ‘This is in contrast with 
the free-streamline model, which leads to a finite drag coefficient (the work 
done by the body appearing as an addition to the infinite amount of kinetic 
energy—using now axes fixed in the fluid at infinity upstream—possessed 
by the fluid). ‘The prediction of finite drag coefficient by the free-stream- 
line model has been greeted as a strong argument in its favour, since it avoids 
the conflict with experience usually described as the d’Alembert paradox*. 
My own view is that steady laminar flow in the limit of infinite Reynolds 
number is so far outside the range of ‘experience’ that the prediction of a 
finite value of the drag coefficient under these conditions is not more wel- 
come physically than the prediction of a zero value. Indeed, in view of the 
indisputable result that the drag coefficient for a body of streamline form, 
for which the development of the laminar boundary layer over the entire 
boundary is calculable, is of order R-!? when R is large, there is a slight 
advantage in the prediction of zero drag coefficient for a bluff body in the 
limit R —> ©, since it reduces the distinction between bodies of different 
shape. 

As mentioned already, an investigation of the character of steady flow in 
a closed region at large Reynolds number has been described in a previous 
paper (Batchelor 1956), and the results of this investigation can now be used 
in a description of the flow in the closed wake bubble for certain kinds of 
bluff body. It was found there that, given that the flow is truly steady 
however large the Reynolds number may be, the vorticity distribution in a 
closed region could be worked out in cases of (a) two-dimensional flow, 
(6) axisymmetric flow without azimuthal swirl, and (c) axisymmetric flow 
with swirl, subject to some restrictions. It is not necessary to describe all 


*Inasmuch as in real flow past a bluff body the drag coefficient is finite—due, in 
all probability, to the existence of turbulence in the wake—the free-streamline model 
doubtless provides a closer representation of the general form of the pressure distri- 
bution on the body in the real flow than does the classical solution which assumes 
wholly irrotational motion. However, the representation of certain aspects of the 
real turbulent flow is an ad hoc objective, different from that pursued here. 
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the results here; suffice it to say that in case (a) the vorticity is uniform 
(w = we, say), and in case (b) the azimuthal component of vorticity (which 
is the only non-zero component) is equal to «r, where r denotes distance 
from the axis of symmetry and x isaconstant. ‘These vorticity distributions 
are valid throughout a closed region, except in the neighbourhood of the 
surrounding singular surface where viscous forces are not small. The 
constants w, and x, which measure the strength of the inviscid ‘standing 
eddy’ in each case, are determined mathematically by the condition that 
steady motion should be possible in the surrounding viscous boundary layers. 

As an example of the application of these results to the proposed limit 
flow with a closed wake bubble, consider the case of a two-dimensional flat 
plate set broadside-on to the stream. The surface of the plate is one of the 
singular curves, in the neighbourhood of which viscous forces are appre- 
ciable however large the Reynolds number may be, and there will be other 
such singular curves in the fluid coincident with the streamlines passing 
downstream from the edges of the plate. ‘These streamlines mark the 
boundary of the wake bubble, which we are supposing to be closed at the 
downstream end, and their shape is shown schematically in figure 1. ‘The 
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Figure 1. Flow past a two-dimensional flat plate in the limit of 
infinite Reynolds number, as proposed herein. 
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motion inside the wake bubble is driven by the action of viscous stresses at 
the curve AC, so that the streamlines just outside AC have greater total head 
than those just inside AC. In the limit R-> ©, the streamline AC will 
thus be a vortex sheet, whose strength will vary along its length according 
to the requirements of continuity of pressure across AC, viz 
(U/U,)?-—(V/U,)* = constant, h say, 

where U and V are the velocities just outside and inside the vortex sheet 
respectively and U, is the speed of the uniform stream. At the closure 
point C, V is necessarily zero, and this relation shows that U is then finite 
at C. ‘This requires the wake bubble to have a cusp at C. 

The property of symmetry about the line bisecting the plate shows that 
the wake bubble must be divided into two halves, the vorticity taking oppo- 
site signs in the two halves. ‘This suggests that the line of symmetry DC is 
a singular curve, as can also be seen to be the case from a consideration of 
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the viscous shear layer in the neighbourhood of the curve AC. The stream- 
line that intersects the plate boundary at A divides the unclosed streamlines 
outside the wake bubble from the closed streamlines inside the bubble, and, 
since viscous forces are significant in a thin layer extending on both sides of 
the curve AC, it follows that at C part of the boundary layer about AC 
passes downstream and part returns to the plate along the line of symmetry 
CD. ‘Thus the wake slit extending downstream from C is a singular curve, 
so too is the portion of the line of symmetry CD, and finally so too is the 
portion DA of the boundary in view of the need to satisfy the no-slip condition 
there. We are led to a simple picture of the wake bubble as a pear-shaped 
domain in which viscous forces are negligible everywhere except near the 
boundary and the line of symmetry, the inviscid flow not near these curves 
being characterized by uniform vorticity, of strength — wy, in the upper half 
of the bubble and +, in the lower*. (A propos of the point made earlier, 
about the possible value of a knowledge of the limit flow as a means to the 
calculation of the flow at finite Reynolds numbers, it will be noticed that the 
general form of the flow past a circular cylinder at Reynolds numbers 
between 20 and 30 (see Goldstein 1938, §20), and that at Reynolds 
number 40 (as found numerically by Kawaguti 1953), both resemble that 
of the corresponding limit flow proposed here.) 

The closed line ACDA isa streamline lying within a continuous boundary 
layer surrounding one half of the wake bubble. Over the part AC of this 
closed curve the viscous stress exerted by the external flow tends to produce 
a clockwise rotation of the upper standing eddy, on the part DA the stress 
opposes circulation in any direction, and on the part CD the stress is zero 
owing to the symmetry. The standing eddy takes up an equilibrium rate 
of rotation, measured by the vorticity —w », in response to these driving 
forces; the larger is the ratio of the length of boundary over which the eddy 
is driven to that over which there is retardation, the larger will be w,, and 
the smaller will beh. An equivalent statement is that w, takes such a value 
that the relative momentum of the boundary layer about ACDA is augmented 
over part of its path and depleted over other parts, the net effect being to 
allow the boundary layer to return to its starting point without change. 

It is not difficult to see how this picture of the flow should be changed to 
suit bodies of different shape. Ifa two-dimensional plate is not broadside-on 
to the stream, the flow has the same general features as for the broadside-on 
position, except for the symmetry; in particular, the dividing free 
viscous layer CD is no longer straight, and the vorticities of opposite sign in 


* Actually, this picture may be a little too simple since it is possible that the 
fluid at the centre of the boundary layer on CD is brought to rest, before reaching D, 
by the adverse pressure gradient near the stagnation point at D, and that the free 
boundary layer cuts across the corner at D. In this event, another singular surface 
would exist to divide the main body of fluid in the wake bubble from a secondary 
standing eddy in the corner at D, and indeed there may even be a whole sequence of 
such singular surfaces and standing eddies of diminishing size as D is approached. 
These are refinements to the picture and I propose to ignore them in order not to 
obscure the main argument. 
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the two sections of the wake bubble are no longer of equal magnitude (but 
note that the loss of head across the bubble boundary, represented by h, 
must be the same for the two sections of the bubble in view of the need for 
continuity of the pressure across the dividing free layer CD at position C 
where the velocity is zero on both sides). In the case of two-dimensional 
bodies without sharp edges, there is the additional complication that the 
two positions of separation of the boundary layer‘on the forward portion of 
the body are no longer determined by the geometry alone, but in other 
respects the picture is unchanged. In the case of bluff axisymmetric bodies 
placed symmetrically in the stream, the singular surface in the interior of 
the wake bubble degenerates to a line, and the free viscous layer returning 
from the closure point along the axis of symmetry in the form of a thin 
cylinder no longer separates two regions of inviscid motion. ‘The distri- 
bution of vorticity within the wake bubble is now given by w = ar, and 
the bubble boundary has the same general form (schematically) as that in 
figure 1. Again the constant x is determined by the condition that the 
motion in the viscous layer bounding the wake bubble should be steady. 
For a three-dimensional body without axial symmetry, the same general 
principles apply, but the form of the distribution of vorticity within the 
wake bubble is not known. 

The proposed model of steady flow at very large Reynolds number has 
now been described, and there remains the question, does such a flow 
satisfy the Navier-Stokes equation? ‘There are two parts to this question, 
one concerning the flow of the thin layers where viscous forces are signi- 
ficant, and one concerning the regions in which the flow is essentially inviscid. 
Consider first the viscous layers. One viscous layer forms on the forward 
portion of the body, separates from it and passes downstream in contact 
with the wake bubble, leaves the wake bubble at the closure point, and 
passes to infinity downstream as a wake. ‘The development of this viscous 
layer is not subject to any restrictions, and will proceed according to what- 
ever inviscid velocity distribution exists just outside the layer. Another 
viscous layer has a closed path, which includes the back of the body and the 
inner side of the boundary of the wake bubble (being in contact with the 
other viscous layer over this part of its path). ‘This viscous layer is subject 
to the restriction imposed by the need to return back on itself after one 
circuit of its path, and steady viscous flow will be possible only if the inviscid 
velocity distribution just outside the layer has a suitable form. I think it is 
reasonably certain that this restriction can be satisfied by an appropriate 
choice of the intensity of the inviscid motion within the wake bubble, the 
latter being measured by the constants w, and « in the simple cases described 
above. If the boundary of the wake bubble did not change its shape with 
change of w, or x, this proposition would be beyond question, and it seems 
likely to remain valid—although no longer readily provable—when the 
changes of shape in the bubble are admitted. 

The flow in the regions of inviscid flow must now be shown to be self- 
consistent, bearing in mind that, for a given shape of the bubble boundary, the 
inviscid flow inside the wake bubble is already prescribed—as to distribution, 
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by the requirements of steady flow in a closed region, as worked out in 
detail for certain simple cases in the earlier paper, and as to intensity, by 
the need for steadiness in the surrounding viscous layer. Moreover, the 
flow outside the combined body and wake bubble is prescribed by the 
irrotationality, as soon as the shape of the bubble boundary is given. ‘Thus 
the crux of the inquiry lies in the possibility of finding a bubble shape, for 
each of a number of arbitrarily chosen values of wy or « (speaking of simple 
symmetrical bodies for the sake of definiteness), and for a given position of 
intersection of the bubble boundary with the body surface, such that the 
inviscid motions in the two sides of the bubble boundary satisfy the condition 
(U/U,)?-(V/U,)? = h, 

where / is an unprescribed constant (which may vary with w, or «). If 
such a family of bubble shapes for different values of wy or « within some 
continuous range can be found, then presumably for one of them the con- 
dition for steadiness of the closed viscous layer will be satisfied and the 
Navier-Stokes equation will be satisfied everywhere. 

So far as the regions of inviscid flow are concerned, the proposed model 
of the limit flow can be regarded as differing from the free-streamline 
model by its allowance for motion inside the wake bubble (as well as by its 
rejection of open bubbles). As stated in § 2, it is generally believed that the 
free-streamline model for a given body is determinate when either the 
position of intersection of the free-streamlines with the body surface or the 
wake pressure is given, the required condition coming presumably from a 
consideration of the boundary layer on the forward portion of the body 
surface. ‘The model proposed here replaces the wake pressure parameter 
by the parameter h (the two being uniquely related when there is no motion 
in the wake bubble), and introduces an additional parameter w, or «, charac- 
terizing the motion inside the wake bubble, which must be determined from 
a consideration of the closed boundary layer surrounding the wake bubble. 
It seems a reasonable inference from experience with the free-streamline 
model that, for a given body shape and position of intersection of the bubble 
boundary with the body surface, the value of h is determined by the choice 
of w) or «. A ‘relaxation’ solution of the inviscid equations for the two- 
dimensional flow down a finite step in an otherwise plane boundary, the 
region in the lee of the step being occupied by a standing eddy with uniform 
vorticity wp», has been obtained by Miss Ann Hawk at Cambridge for one 
arbitrarily chosen value of wy, and the value of # was here found to emerge 
as one of the results of the calculation. 

Problems of inviscid flow with finite vorticity and boundaries of un- 
known shape are very difficult to handle, and I have not been able to show 
in general that it is in fact possible to find a bubble shape for a given body 
shape and position of the intersection of the wake bubble with the body 
surface, and an arbitrary value of w, or. A little support for the belief that 
an appropriate bubble shape exists for each value of wy or « (within a certain 
range) can be found in the fact that a Hill’s spherical vortex (see Lamb 1932, 
ch. 7) has a vorticity distribution of the required form, viz. w = ar, and 
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satisfies the above boundary condition with h = 0 when « = 15U,/4a? 
(a being the radius of the sphere containing the vorticity). A Hill’s spherical 
vortex may be regarded as a possible wake bubble for a body in the form of 
a sector of a spherical shell subtending an angle 4), say, at the centre of the 
sphere, so far as conditions outside the viscous layers are concerned, « 
having the maximum possible value and h being zero. At the other extreme 
there is the free-streamline solution for the same spherical cap (the details 
of which are not known mathematically, although such a flow has been 
observed about a cavitating sphere in water), corresponding to « = 0 and to 
whatever is the value of 4 that corresponds to the intersection of the free- 
streamlines and the body being at the edge of the spherical cap. For all 
intermediate values of x, and consequently of f, there should exist an in- 
viscid wake bubble for the same spherical cap, the bubble boundary having a 
finite length (at any rate, for « 40) and a cusp at the rearmost point. ‘The 
value of x for which the viscous layer bounding the wake bubble is in steady 
motion seems likely to be one corresponding to a bubble whose length is 
not many times larger than the linear dimensions of the spherical cap—for 
if it were much larger, the area of bubble boundary over which the irro- 
tational stream accelerates the standing eddy would be much larger than 
the area over which the cap retards it, and the vorticity in the standing eddy 
would take up a value close to its maximum, and this in turn corresponds 
to a short bubble. 
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SUMMARY 
The compression of a perfect gas between a uniformly moving 
piston and a rigid wall is discussed in the one-dimensional case. 

If the piston moves with a finite speed, it will initiate a shock in the 

gas which will reflect successively from rigid wall and piston 

and cause the compression process to deviate from a reversible 

adiabatic process. Expressions are derived for the relative 

changes in pressure and density at each shock reflection. Then 
values of density and pressure after any number of shock reflections 

are computed relative to their initial values, and, in terms of these, 

the corresponding values of temperature and entropy, as well as 

shock speeds, are determined. ‘The limiting value of the entropy 

change, as the number of reflections goes to infinity, is obtained as 

a function of the ratio of specific heats of the gas and the strength of 

the initial shock. Hence it is possible to estimate an upper limit 

to the deviation of the shock compression process from a reversible 

adiabatic process. Some illustrative numerical examples are 

given. 
I. INTRODUCTION 

When a gas is compressed by the motion of a piston, one phenomenon 
that may cause the compression to deviate from a reversible adiabatic 
process results from the fact that, if the piston moves with a finite speed, 
hydrodynamic shocks will be produced. It is the purpose of this paper 
to discuss, in the idealized case of zero heat conduction and viscosity, 
the compression of a gas by a piston which moves with finite speed, to give 
a detailed description of the process as shocks run back and forth between 
the moving piston and an enclosing rigid wall, and to determine the deviation 
of this shock compression process from a reversible adiabatic process. 
In the absence of dissipative effects, the shock front may be represented 
as a mathematical discontinuity in pressure, density, etc.; for a discussion 
of this point, see Courant & Friedrichs (1948), also Grad (1952) and Sachs 
(1946). 

Consider the one-dimensional system consisting of a perfect gas with 
zero viscosity and heat conductivity confined between two plane parallel 
walls that are impermeable to heat. One wall is considered to be rigid 
and the other to move toward the first with constant speed. If the speed 
of the moving wall (or piston) is infinitely slow, then (assuming the internal 
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energy of the gas to be proportional to the temperature) the pressure p 
and the density p will increase continuously according to the adiabatic 
relation pp” = constant, where y is the ratio of specific heats, and the 
process will be reversible. If, however, the piston is suddenly given a 
finite velocity inward, it will initiate a shock in the gas. The shock will 
reflect at the rigid wall producing a second shock which will in turn be 
reflected at the moving wall, and so on. When each new shock passes 
a given point, the pressure and density of the gas at that point will undergo 
discontinuous increases. As a shock front traverses the distance between 
piston and rigid wall, it separates the gas into two regions throughout each 
of which the pressure, density and material velocity are uniform, but change 
discontinuously across the shock front. At the instant of reflection at 
either end, one of the regions has zero volume and the gas is uniform 
throughout. Since the entropy of the gas increases as each shock passes, 
the process is irreversible and pp~’, which is a function of the entropy, no 
longer remains constant. As a consequence, at a given compression the 
pressure and temperature of the gas will be higher in the shock compression 
case than they would be in the case of constant entropy. ‘The amount 
by which the process deviates from the adiabatic relation pp~’ = constant 
depends upon the strength of the initial shock ; that is, the ratio of the pressure 
behind the initial shock to its original value, which in turn is related to the 
speed of the moving wall. 

Using the conservation equations for shocks, it is possible to compute 
the ratios of density and pressure behind the shock to their values in front 
of it after each successive reflection. ‘Then the values of density and pressure 
behind the shock after any number of reflections are computed relative to 
their initial values, and, in terms of these, the corresponding values of 
temperature and entropy, as well as the shock speeds, are determined. 
In order to estimate an upper limit to the deviation from the constant entropy 
case, the limiting value, as the number of reflections goes to infinity, of the 
function pp” is computed. ‘This limit is derived as a function of the strength 
of the original shock and the y which characterizes the gas. Finally, some 
illustrative numerical examples are given. 


I]. PRESSURE RATIO AND COMPRESSION PRODUCED BY EACH SHOCK 

Let p,, p, and u,, (nm = 0,1, 2,...) be the density, pressure and material 
velocity, respectively, behind the mth shock, where the condition n= 0 
corresponds to the initial conditions in the gas, m = 1 corresponds to condi- 
tions behind the first shock (that is, the initial shock produced by the moving 
wall), m = 2 corresponds to conditions behind the second shock (that is, 
the shock produced by the first reflection), and so on. The material 
velocities, “,,, will be measured relative to the rigid wall and the gas is assumed 
to be initially at rest; that is, v4, = 0. Let 


On=PnlPn—w» (1 a) 


Nn = Pn|Pn-1 (1b) 
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The equations describing conservation of mass and momentum may be 
combined to give the relation* 
ae —— 
(U,4 Ge) (Pn —Pn- (Pr-1 —Pn ). 
The initial and boundary conditions are such that 
(a) if the shock is produced at the rigid wall, 
u, = 0, U,41 = W; 
(b) if the shock is produced at the moving wall, 
U, = UW, u,-1 = 9, 
where w is the velocity of the moving wall. So in either case the left-hand 
side of the above equation is w*, and, substituting (1), it becomes 
gz x a —n-! 
W"Pn—1 Pri = (c,, 1)(1 Qn ). (2) 
If the perfect gas has the property that the internal energy is proportional 
to the temperature, then the internal energy per unit mass is given by 
E, = p,lp,(y—1), 


and the conservation equations lead to the Hugoniot relation 





(y—1)-—(y+1)y, 1—p7,, 
on = } : ( a ’ 3 
(~~ Oth nF \ ) 
where B= (y+1)/(y-1), 


and y is the ratio of specific heats. 

To obtain expressions for o,, and 7,, as functions of m, y and oj, first 
write (2) with m replaced by n+ 1, and divide the resulting equation by (2). 
The result may be written 

‘ + -1 

Go lh +1 1)(1 — n+ 1) on Ni(Fn- 1)(1 — An ). 
Then, by substituting (3) into this equation, one obtains the equation 

(u Vn 1)(y n+1— 1)? = (n _ 1)°(u— Nn+1)In+1 
which is a quadratic having the two solutions: 

_ y-1. 
Qn+1 =n > 

and Anat = (HIn—1)/(nn +42). (4) 
The first of these solutions does not apply because the compression, 7,,; 
must be greater than or equal to unity for all 7. 

Given the strength of the initial shock (that is, o,) and using equation (3) 
together with the recursion formula (4), it is possible to obtain o, and 7, 
for all. ‘The general expressions are found to be 

o, = (A+ptn)/(A+n-—1), (5a) 
n, =(Atptn—1)/(A+n), (5b) 

where 
A=(uH+1)/(,-1), w=(y+)/(y-1), 1<y<a, 1<a<o. (6) 

* For a derivation of the shock relations and a general discussion of shock waves 
in one-dimensional flow see, for example, Courant & Friedrichs (1948), Chapter ITI, 
Part ©. 


F.M. 2D 
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That equations (5) do indeed represent the general solution may be verified 
by induction as follows: 
(1) Form = 1, (5b) becomes 7,, = (1+ 20,)/(u4+,), which, on comparison 
with (3), is seen to be the correct result. 
(2) Equation (5 b) satisfies the recursion relation (4). 
(3) Equations (5) together satisfy (3). 


III. CONDITIONS BEHIND THE 2TH SHOCK 


1. Pressure and density 
Let z,, and «,, be the pressure and density, respectively, behind the mth 
shock referred to the initial pressure and density, p) and py. ‘Then, using 
(1) and (5), one obtains 
(A+pt 1j(At+p4+2)...(At+ut+n) 








7™n = PniPo — _ AA + 1)...(A+n— 1) , 9) 
(A+p)(A+p+1)...(Atp+n—1) - 
ae = 
Kn = PniPo = bide. (A+ 1)(A+2)...(A+n) a 
On introducing the difference equation for the ['-function, namely, 
I(z+1) = 2I(z), (8) 


(7) may be written 








_ PATA+H+n41) 

7a = TAt+p+llatn)’ os 
— PA+1)PA+pu+n) 

“o™ TA+plata+i)’ mm 


2. Temperature 
Let 0, = p/p. Then the temperature behind the mth shock referred 
to its initial value is, from (7), 
6, - Ty A+n A4 stn 
ee ee ee ae ol a (10) 
oh Po Pn Kn NA +) 





3. Entropy 

Rather than deal with the entropy directly, consider the following 
function of the entropy (from which, if desired, the entropy can readily 
be obtained): FF (S) = pp-” = exp[(S—S;)/C,] 
where S is the specific entropy, Sj is a constant and C, is the specific heat 
at constant volume which is assumed to be constant. Let F, = F(S;). 
Then the value of this function of the entropy after 7 reflections, relative 
to its initial value, is given by 


€, = F,,/ Fo = exp[(S, — So)/Cy] a (Pn/Po)(PolPn)” 


jae n 
7 C 








vs = 
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The substitution of (10) and (9b) into (11) leads to 
_ tn _ Atm A+ut+m)(TA+MEA+a+)P4 A, 
oa (A+ pz) PA+1)PA+p+n) : 


The limiting value of €,, as n—> o will be considered in Section IV. 








4. Shock speed and time 
Let ¢,, be the time at which the nth reflection occurs, that is, the time 
at which the (n+ 1)th shock is produced; and let tj =0. Then at time 
t,, the gas has a uniform density equal to p,, and the compression relative 
to the initial density is p,,/py =«,. Let L be the initial distance between 
the piston and rigid wall. Then the distance between piston and rigid 
wall at time t,, is L—wt,, where w is the speed of the piston, and the com- 
pression at time ¢,, is given by 
L 
“a= Tat,’ 
If T = L/w is the time necessary for the piston to traverse the distance L, 
the above equation yields 
t, = T(1-«,?). (13) 
The time interval between the mth and the (m+1)th reflections, that is, 
the time between the (m + 1)th and (m+ 2)th shocks, is given by 
Atnsie = ti -t= T(k,* — Ka) 
From (7b), (1b) and (5b), it follows that 
Knit = KnPn+ 1/Pn = Kia = K,(A+ p+ n)/(A+n+ 1); 
and hence 
Atniue = Tx, (u—1)/(A+p+n), (14) 
where x, = 1. 
Let D,, be the speed of propagation of the mth shock. If m is odd, the 
nth shock originates at the moving wall, and 


_ L-ut,_, 
n(odd) ~ a ay “ 
On the substitution of (13), (14) and L = wT, this equation becomes 
Dyyoaay = (A+ n)/(H— 1) + 1]. 
If n is even, the nth shock orginates at the rigid wall, and 
L—ut,_,—waAt,_+j2 


) | ae = At ; = 
n—1'2 





In the same way as above, this equation becomes 
iitiman a wA+n)/(u— 1). 
The alternative expressions may be combined in the form 


‘ wer if n is even, (15) 
n = WwW +a , = . ° 
oes ) 1, if mis odd. 





2D2 
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IV. LIMITING VALUE OF THE ENTROPY CHANGE 
Let €°, be the limiting value of e, as the number of reflections becomes 
infinite. ‘Then, from (11), 


, : se . ae 
€x (y, 9) = lime, = lim — = II =, (16) 
n> x n>o Fy i=1 i 


and, by the substitution of (12), 
1 R(A+ pu) ]-1 ( P(A+n+1)]') 
0 ' : } ! I t . 
e? (y,o,) = ———-| ——— lim <(A+n)(A+p+n)| ————— P 
(Ys 01) UsALTO DH] sm Crmetet posure) f 
3y the introduction of the asymptotic formula for the gamma function 
T(z) ~ (27)¥2e-727-12, (17) 
and by the use of (6) and (8), the limit on the right of the above expression 


becomes 


e’+1 Jim [(A+n)/(A+p+n)]?-™ = 1, 


noe 
where the final limit is evaluated by taking the logarithm and applying 
l’Hospital’s rule. Hence, 


0: £. a. 1 om 1A T [t) 1 
Cy (y, 7%) a NA \ min \ > | , (18) 


This expression may be generalized to give the limiting value of F(S) 
measured relative to its value behind any shock, say the Ath, where is 
finite. Thus let 








ny De 5G, 
e* (y,0,) = lim c= es, ke = 0, 1, 2,.... (19) 
n>a k t=k+1 ‘ii 
Then, using (5), one obtains the recursion formula 
k wn  O, os (A+k-1)A+p+k-1\" as 
oe = lI - = — € ity = - cy, ’ 
Op ¢=k 1; a ee: (A- pet R)(A : ky 


from which, together with (18), it is readily demonstrated by induction 


that 
: l P(A+u+k))7-1 
fe - 
colo) = ATRaAT es 5 Tas kh > | (20) 


It is of interest to obtain the limiting values of &, (y, o,) corresponding 
to the extreme values of y. Thus, as y > 1, then p > 0, A—> o; and, 
because of (17) and (6), 

e (1,0,) = lim &,(y,0,) = lim e 2D (A+ w+k)/(A+ Ry] OA+2 etd 


y>1 > x 





wok ji. lalla (21) 
which is seen to be independent of k. In the case where y — ©, then 
p> 1, AA) = 2/(¢,-—1); and by taking the logarithm of the second factor 
in (20) and expanding it into a power series in (u—1), one obtains 


aes ! 21"(\y +k +1) 
ee Pes a = 2X an 
Fol ay = Rel a eG yes PA) +k+1) 3 ( 


v= o A) 





7 


(Tables of the function ¢(z) = I'’(z)/P'(z) are given in Davis (1933).) 
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In the limiting case where the initial shock is a “strong shock ’* (that is, 
where o, —- © and 7, — pz), it is seen from (6) thatA—+0. Hence, the 
preceding formulae may be reduced to those applying to the special case 
where the initial shock is a ‘strong shock’ by setting A = 0. It is seen that 
(20) becomes infinite for this case if either k = 0 or y = 1, but otherwise 
remains finite. 

V. NUMERICAL EXAMPLES 

In order to compare the case of compression at constant entropy with 
the process of shock compression discussed in the foregoing sections, let 
the same symbols as used above, but without subscripts, represent corres- 
ponding quantities in the case of constant entropy. ‘Then the pressure z, 
the density « and the temperature 7, referred to initial values, are related by 
the expressions 


1 
TK IF 7K 1 = 7, 


The corresponding relations for conditions behind the mth shock are 


Tak. = 'e,, Wig Key Tig 


Thus, with the same initial conditions and the same compression x = k,, 
the relation between the pressures that would be arrived at by the two 
processes is 

PrlP = Trl = En; 


and, since temperature is proportional to pressure at a fixed density, 


As the number of reflections becomes large, €,, approaches the limiting 


n 
x, are confined toa strip in the (z,«~!) plane (PV diagram) which is bounded 
by the two adiabatics z«~” = 1 and a«~” = e€, (y,o,).. The case of a diatomic 
gas (y = 1-4) is illustrated in figure 1 for two different values of the strength 
of the initial shock, o, = 5ando, = 50. The figure representsa PV diagram 
and «~! are pressure and specific volume 


value €”, (y,o,). Hence, for a given y and o,, the possible values of 7, and 


=] 


on a logarithmic scale, where 
expressed in units of the initial pressure and specific volume, py and pj’, 
respectively. The points (7,,*,,), plotted as circled dots, are seen to be 
converging quite rapidly to the adiabatics for the limiting values of the 
entropy. In table 1 numerical values of €}, (y,o,) are listed for several 
different values of y and o,. 

For a given y and oj, the greatest relative change in F(S) = pp~” caused 
by the various shocks occurs in the first shock. It is of interest to see to 
what extent the process of shock compression subsequent to the first 
shock will deviate from a constant-entropy process that starts from condi- 
tions corresponding to those produced by the first shock. In the limiting 
case where the first shock is a ‘ strong shock’ (o, = ©), the first shock causes 

* The special case where the initial shock is a ‘ strong shock’ has been discussed 


previously in an unpublished report by the present authors (1953). 
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an infinite change in entropy, but the subsequent change as the number 
of reflections goes to infinity remains finite for y > 1. In table 2 numerical 
values of e', (y,o,) are listed for several different values of y and o,. In 
figure 2a PV diagram similar to that of figure 1 is shown, where now z and 
x! are pressure and specific volume expressed in units of the pressure and 
specific volume behind the first shock, p, and p;', respectively. The case 
illustrated is that of a ‘strong shock’, which, for a given y, gives the greatest 
deviation from constant entropy. 
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Figure 1. PV diagram for y = 1-4, where 7 and «~! are pressure and specific 


volume in units of their initial values. The straight lines are adiabatics for 
the initial value of entropy (solid) and the limiting values of entropy (dashed) 
in the two cases: o, = 5; o, = 50. The circled dots represent values of 
pressure and specific volume behind the mth shock referred to their initial 
values, where the numbers indicate values of m in each case. 
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\ “ae 1-4 5/3 2 3 00 
aN 
2 | 1-083 | 1-077 | 1-074 | 1-071 | 1-067 | 1-055 
5 | 1513 | 1511 | 1:506 | 1-500 | 1-485 | 1-434 
10 | 2:26 2-31 2-32 2-31 2-28 2-18 
25 4-42 4-72 4-78 4°81 4:76 4°52 
50 | 7-94 8-75 8-88 8-94 8-91 8-45 
100 | 14:85 | 16-83 | 17-14 | 17-27 | 17-25 | 16-33 





























Table 1. a as a function of y the ratio of specific heats of the gas, and o,, the 


pressure ratio in the original shock. €®, is the value of the ratio of pp~” after 
an infinite number of reflections to its initial value, where p and p are pressure 
and density, respectively. 
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Figure 2. PV diagram for y = 1:4, where 7 and x~! are pressure and specific 
volume in units of their values behind the first shock. The straight lines are 
adiabatics for the value of entropy behind first shock (solid) and limiting value 
of entropy for the case o, = ©, i.e. a ‘ strong shock’. The circled dots 
represent values of pressure and specific volume behind the th shock referred 
to their values behind the first shock, where the numbers indicate values of n. 
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. } 
\ : 5/2 ) 2 - 
ae 1 1-4 5/3 2 3 | S 
Z 1-083 | 1-062 1056 | 1-050 | 1-042 1-028 ) 
5 TPS 1-29 1:24 | 1:20 | 1:15 1-09 
10 | 2-26 1-50 1:39 | 1:31 | 1:22 | 1-12 
25 | 4-42 1-73 1°53 142 | 1:28 | 1:14 | 
50 | 7-94 1-84 1-59 1-45 | 33 | 1-15 
100 | 14-85 1-91 1-63 148 | 1:32 | 146 | 
s / 1-9 1-66 1:50 | 1:33 | 1:16 | 
Table 2. «1, as a function of y, the ratio of specific heats of the gas, and o,, the 
pi ratio in the original shock. 1, is the value of the ratio of pp~? after 
ite number of refk ms to its value behind the first shock, where 
ire pressure and sity, respectively 
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The method of characteristics for steady supersonic 
rotational flow in three dimensions 


By MAURICE HOLT 
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SUMMARY 

The method of characteristics for steady supersonic flow 
problems in three dimensions, due to Coburn & Dolph (1949), is 
extended so that flow with shocks and entropy changes may be 
treated. Equations of motion based on Coburn & Dolph’s 
characteristic coordinate system are derived and a scheme is 
described for solving these by finite differences. 

A linearized method of characteristics is developed for cal- 
culating perturbations of a given three-dimensional field of flow. 
This is a generalization of the method evolved by Ferri (1952) for 


perturbations of plane flow and conical flow. 


1. INTRODUCTION 

The method of characteristics for problems of steady supersonic flow in 
two independent variables has a sound and clear cut basis which appears to 
be lacking when additional independent variables are considered. In two 
dimensions two distinct families of characteristic curves are defined uniquely 
and the independent variables can be transformed so that the problem 
expressed in terms of characteristic coordinates is equivalent to the original 
problem. ‘The extension to fields involving three dependent variables 
presents no difficulty. ‘Vheorems of uniqueness and existence covering 
most types of two-dimensional flow have been established. 

Flow involving three independent variables has a more nebulous 
character. A one-parameter family of characteristic surfaces pass through 
any one point in space and these envelop a conoid through the point; the 
lines of contact of the surfaces and the conoid, the bicharacteristic curves, 
also form a one-parameter family. ‘There is, therefore, at first sight, a 
certain arbitrariness in any numerical method based on characteristics 
enabling one to determine the flow off a given initial surface. 

Three distinct methods for three-variable problems have been proposed 
so far. ‘Thornhill (1952) proposes two difference schemes, both starting 
from a triangle drawn on an initial surface. In the first a new point in the 
flow space is determined as the intersection of the three internal characteristic 
surfaces through the sides of the triangle. Inthe second the common inter- 
section of the three internal bicharacteristics through the vertices is found. 

Sauer (1950) reduces the three-variable problem to a series of two-variable 
problems. Working in Cartesian space he selects a sequence of equally 





410 M. Holt 


spaced coordinate planes, x = constant, say. In each plane in turn he uses 
the two-dimensional method of characteristics to calculate the network inter- 
sected by families of characteristic surfaces passing through curves 
/ (vy = constant) drawn in the initial surface. Each point in this network is 
connected with a corresponding point in the network in an adjacent plane by 
a difference equation in the x-direction. 

Coburn & Dolph (1949), in a paper of considerable importance, consider 
three-variable problems from a more formal standpoint. They draw 
attention to work by Titt (1939) on general non-linear hyperbolic equations in 
three independent variables. ‘Titt introduces a characteristic coordinate 
system in terms of which the original initial value problem can be reduced toa 
two-dimensional problem and so establishes uniqueness and existence 
theorems for correctly posed three-variable problems. Coburn & Dolph 
suggest that any difference scheme for supersonic flow problems should be 
closely linked with ‘Titt’s coordinate scheme. ‘The present author is in 
agreement with this view. By adopting the Titt scheme, one is sure of 
satisfying the requirements to ensure uniqueness and existence. At present 
no other scheme of equivalent soundness is known. 

Coburn & Dolph have developed these ideas for the equations of steady 
supersonic, homentropic flow, and their characteristic coordinate system is 
defined as follows. A family of non-intersecting, space-like curves are 
drawn on the surface bearing the initial Cauchy data, and the characteristic 
coordinate system is based on the two families of characteristic surfaces 
passing through these initial curves and on the family of surfaces determined 
by the corresponding bicharacteristics. These three families of surfaces are 
taken as coordinate surfaces, and the equation of potential flow is replaced by 
two characteristic equations and symmetry conditions. ‘The characteristic 
system of equations is in a form suitable for solution by a method of finite 
differences but Coburn & Dolph do not consider the details of this. 

Of the methods proposed for practical application Sauer’s is the closest 
to the formal approach of Coburn & Dolph, since he employs two families of 
characteristic surfaces and a third non-characteristic family. However, 
Sauer chooses the third family to be planes x = constant, and in so doing is 
unable to satisfy Coburn & Dolph’s requirement that two coordinate 
directions shall be bicharacteristic. When this condition is satisfied the 
third family of coordinate surfaces is determined automatically. When 
Coburn & Dolph’s method is applied, in general, none of the coordinate 
curves in the characteristic system are known in advance and they have to be 
determined, step by step, together with the physical variables arising in the 
problem in question. 

In the present paper Coburn & Dolph’s work is carried a stage further. 
Firstly, their method is generalized to take account of motion with entropy 
changes and vorticity (viscosity and heat conduction are neglected). Even 
the simplest three-variable problems involve shocks of variable curvature, so 
this step is essential. The basic characteristic coordinate system is set up 
exactly in the manner proposed by Coburn & Dolph. When the motion is 
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referred to this system two of the transformed equations are generalizations 
of the second characteristic equations familiar in non-homentropic problems 
in two variables. A third equation relates derivatives of pressure and velo- 
city components. The fourth equation, the energy equation, directly 
connects total speed with variables of state. ‘The entropy does not enter 
these equations explicitly, and appears only in the equation of state and the 
condition that entropy remains constant in the stream direction. 

A finite difference scheme is considered for the numerical solution of the 
equations in practical cases. ‘The basis of the scheme is the construction of a 
sequence of linear characteristic networks in the surfaces determined by the 
bicharacteristics. ‘The method of construction of any one network is similar 
to that adopted in non-homentropic flow in two dimensions. An extra 
difference equation connects any two adjacent networks. 

In the general case, when derivatives with respect to all space variables 
occur non-linearly, the calculation of the coordinate system is completely 
interwoven with that of the dependent variables themselves, and the computa- 
tion required is indeed formidable. However, if we know the solution to a 
certain non-linear flow problem, the equations governing a field of flow 
which is a linear perturbation of this are linear, and can be expressed in terms 
of a characteristic coordinate system defined by the basic flow. Such 
perturbation problems are of wide practical interest and have been examined 
previously, chiefly by Ferri (1952), but also by Ferrari (1936) and Guderley 
(1947). 

Ferri has developed a Linearized Method of Characteristics to apply to a 
number of three-variable problems of this simplified type. He considers 
three-dimensional linear perturbations of known two-dimensional fields of 
flow. Among the applications he considers are axially symmetrical flow 
past bodies of revolution which differ slightly from conical shape, flow past 
bodies of revolution (not necessarily thin) at small angles of yaw, and fields 
of flow in which there is a slight departure from plane flow conditions. In 
all cases his perturbation equations are referred to the characteristic network 
of the basic flow. In the case of plane flow the third independent variable 
is taken to be the Cartesian coordinate normal to the flow plane. In the 
axially symmetrical case, the angle of rotation of a basic meridian plane is 
used. 

Essentially, in Ferri’s Linearized Method of Characteristics, the departure 
point is some appropriate two-dimensional method. Here we shall arrive 
at a Linearized Method of Characteristics from a different approach. We 
shall first establish the equation governing the perturbations of a basic 
three-dimensional flow referred to Coburn & Dolph’s system of characteristic 
coordinates. We shall then consider how these equations simplify when the 
basic flow degenerates into a plane or axially symmetrical flow. In this way 
it is easier to make sure that Coburn & Dolph’s system of equations is 
retained in the form accepted for fully three-dimensional flow. 

It is difficult to compare the equations obtained here with Ferri’s equations, 
since not only are the two approaches different but, also, the dependent 
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variables are handled in a different manner. However, in the special cases 
treated the characteristic coordinate systems are the same and it should be 
possible to deduce one form of equations from the other. 


2. ‘THE EQUATIONS OF ROTATIONAL MOTION IN A GENERALIZED COORDINATI 
SYSTEM 

Let the position of a point in a field of steady, supersonic rotational flow 

be defined by the coordinate vector x‘ with metric tensor g,,, and let u' be the 

velocity vector at this point. Denote the pressure, density and entropy 

p,p and S respectively. We shall neglect viscosity, heat conduction and 


by ; 
radiation. 
Referred to this coordinate system the Eulerian equations of motion are 
l ») 
ft { 
wu, ,+ —p 0, (2.1) 
pP 


together with the equation of continuity, 


u 
guj, + py =0 (2.2) 
p 
In steady flow there is no variation in entropy in the stream direction, so 
that 
ws (0). (2.3) 


lhe equation of state may be written 


p = p(p, S). (2.4) 


We write c (2). (2.6) 


Then, from (2.3) and (2.5), we obtain the relation 


up, = cup ,. (2.7) 
Using this, (2.2) may be written 
th , 
U; i 3P.j 0). (2.5) 
pe 
In supersonic flow real characteristic surfaces are associated with (2.1), 
(2.8), (2.3) and (2.4). Through each point x! pass a single-parameter family 
of surfaces which are characteristic for the velocity vector, the pressure and 
the density. In addition, each surface through x' which contains the stream 
tic for the entropy. 








direction is characteri 


Consider first (2.1) and (2.8). Then the surface , 


s(x') = constant, 


is characteristic if, when we specify on & values of u', p and p, equations (2.1) 


and (2.8) admit solutions with arbitrary values of derivatives of these 
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quantities leading out of the surface. An elementary calculation shows that 
if A‘ is the unit vector normal to = at x’, then the surface is characteristic if 

(AP = c. (2.9) 
This means that there are 00! characteristic surfaces through A‘, which are 
tangential to a cone, the characteristic cone, with axis along the velocity 
direction. ‘The component of velocity normal to a characteristic surface is 
equal to the local speed of sound. We shall be interested only in the nappe 
of the cone downstream of x’, on which 

hear. . 
wr, = +¢. (2.10) 
We observe that these characteristic surfaces are defined exactly as in 
homentropic flow, so we may make use of the properties of such surfaces 
deduced by Coburn & Dolph. 
Now consider (2.3) and write 
S;=Ajst+pjtt+v,;u, 
where X, is a unit vector normal to the surface and ,, v; are unit vectors lying 
in the surface 
2(x') = constant. 

Then we find 

! j er ae od 

wA,;stu’p,tt+u'v,u = 0, 
and we are therefore unable to determine s uniquely if 
wX; = 0. 

It follows that surfaces containing the velocity vector are characteristic for 
the entropy. 


3. EQUATIONS OF MOTION IN CHARACTERISTIC FORM 

To obtain the characteristic equations corresponding to (2.1) and (2.8) we 
may employ the same coordinate system as that defined by Coburn & Dolph 
for homentropic flow. 

In Coburn & Dolph (1949), initial data (in the general case, values of 
u', p, p, S compatible with the equation of state) are given on some non- 
characteristic initial surface. On this surface an co!-family of curves are 
drawn which are space-like with respect to the local Mach cone. ‘Through 
each curve of this family pass two surfaces which are characteristic in 
relation to (2.1) and (2.8). In this way two families of characteristic surfaces 
are defined in the region of flow and these will intersect in an 0*-family of 
curves, which will contain the original co1-family on the initial surface. 

At each point on a curve of this «?-family two bicharacteristic directions 
are defined. These are the directions of the lines of contact between the 
two characteristic surfaces and the characteristic conoid at the point. ‘Thus 
at each point in the flow region three directions are defined. On these 
Coburn & Dolph base their characteristic coordinate system. 

Let /‘ be the unit vector along a line of the oo?-family. Let A’, A“ be unit 
vectors at this point directed along normals to the two characteristic surfaces 
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through the point, and let ¢', t’' be the corresponding bicharacteristic 
directions at this point. The required coordinate system is to be based on 
the triad //, t, t’, 

Following Coburn & Dolph we now transform (2.1) and (2.8) so that they 
are referred to the coordinate triad /', t', t. In particular, we require that 
only derivatives in the new coordinate directions shall appear in the equations 
of motion. 

We express u, ; as the sum of components along X’, t’, J, Then 

u;; = Aya, +t, b+]; ¢, 
and wu, ; = wr,a,+wt,b,+u'l; ¢,. 


Equation (2.1) becomes 

wd, a;+u’t; b+ ul; ¢,+ +i, = 0, (3.1) 

Equation (2.8) becomes ° 
BA, a; t+t,b;+l,¢;) + =a"'P =: O, 


4 


, _ 
or Na; +tb,+Ve;+ oat up, = 0. 3.2) 
‘Take the scalar product of A' with (3.1). Then 
; a ... 
(u’A;)A‘a, + wb, N'b, + wl, Ac, + - A‘p , = 0. (3.3) 
pb 


Multiply (3.2) by ¢ and subtract from (3.3). Then, using the property 


A‘u; = +c (A‘is the outward normal on the downstream Mach cone), we find 


(ut; A'— ct')b, + (w’l; X‘—cl')c, + = (ex u')p , = 0. (3.4) 
If q is the total speed we may write (see Coburn & Dolph, equation (2.6)) 
ul = cr + 4/(q? —c?)t", (3.5) 
and (3.4) reduces to 
(ut, A'—ct')b, + (wl; A —cl')c, — a fp,= 0. (3.6) 
pe ; 


We now evaluate b, andc,. We find 
b+tlc, = bu; ;, 
Pt,b,+6, = Pu, 5, 
so if lt, = cos¢ =a, (3.7) 
b(1—cos*s) = tu, ,—al’u, ;, 
¢(1—cos*s) = Pu, ;—at’u, ;, 


and (3.6) reduces to 





~~" {u’(t; — al,) A — c(t'—al')} t*u,, + 





+ {ui(l, — at,)—c(li—at')} I*u, ,— oy =0. (3.8) 
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Similarly, starting from the resolution 
u,,= Aja, t+tybi+h¢, 


we can derive the equation 





1 , 
ize {u?(t — al;)A"'— c(t" —al')} t’*u, + 





+ = 
1—a? 


fyi(] i c(Li— at’) | VIP = 6) ix 0. (3.9 
{w(l; — at))A" — e(l'—at")} l*u, .- ee ‘Pp, =9. (3.9) 

Equations (3.8) and (3.9) are the two ‘second characteristic conditions’. 
They are generalizations of Coburn & Dolph’s equations and express the 
conditions that derivatives of u' or p (or p) leading out of the characteristic 
surface exist, even though they are not determined uniquely. Equation (3.8) 
contains only derivatives along ¢' and /', and equation (3.9) contains only 
derivatives along ¢’' and /'. 

We require three further equations to determine u', p, p, S in the 
characteristic coordinate system. ‘The first is a further relation between 
derivatives in the ¢', t’' and /‘ directions. 

Equation (2.18) of Coburn & Dolph (1949) applies here and is 
-: : (t'+2')- i 

(1—d) V(q* —c?) (1—a*)V/(q?—c’) 
where d = cos¢ = ff. (3.11) 
Substitute the expression (3.10) for u’ in (2.1), but leave u,; unchanged, 
and take the scalar product of /' with the resulting equation. In this way we 
find 








li (3.10) 


ui 


2 
C2 


(1—d)y (q?—c’) 





li (tu, ,+ tu, ;)— 





7 ac? - 1) _ 
GaaVe a" u; + : 1p, =0. (3.12) 

In most practical applications the flow originates from a region of constant 
total energy and the relation 


gt+ {2 = constant (3.13) 


is true throughout the fluid. The constant on the right has the same value 
throughout the flow field irrespective of the presence of shock waves. 
Equation (3.13) provides the second relation required. The third rela- 
tion involves the entropy. 
Suppose that v‘ is the unit vector normal to the surface determined by / 
and the velocity vector u'’. Let m' be the direction intersected by this surface 
on the surface determined by the bicharacteristics t‘ and t’. Write 


S,=vA+m,B+1,C. 
Then B+mil,C = m'S,, 
I'm,B+C=lS,, 
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and if 2 = I'm, 
1 ‘ 
B= —,(m'S,—el'S,), 
1—eé? 
= 4 
C= —,(l'S,-em'S,). 
1—¢2' : yt) 


Hence, from (2.3), 
wm; B+u'l;C = 0, 


or (u’m,; —eu’l;)m'S ,+(u’l; —ew’m, US , = 0. (3.14) 


The jinite difference scheme 

The equations of motion referred to characteristic coordinates, (3.8), 
(3.9), (3.12), (3.13) and (3.14), are in the form most suitable for solution by 
numerical methods and we now put forward a finite difference scheme. The 
basis of this is the construction of a sequence of linear characteristic networks 
in the surfaces determined by the bicharacteristics. ‘Three difference 
equations are employed in each surface and an additional difference equation 
connects the point to be determined with the corresponding point in the 
adjacent surface, where conditions are already known. 








Figure 1. The cell of the characteristic network. 


The cell in the network is shown in figure 1. We start from the quad- 
rangle 1234 on a space-like initial surface. ‘The segments 23 and 14 lie 
on initial lines (/') on the surface. The segments 12 and 43 lie on surfaces 
determined by bicharacteristics on the initial surface. 

It is assumed that the characteristic network in the surface through 12 is 
known. (The first surface found may be ina plane of symmetry.) 

Segments 75 and 86 are the intersections of the surface through /' and u' 
on the elements 346, 215. 
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Equation (3.8) is applied along the segment 46. The values of t*u; ; 
are taken to be those at point 4, and t“u;,, t*p,, are replaced by difference 
ratios Au, Ap 

distance 46’ _— distance 46° 

Equation (3.9) is used in a similar manner to set up a difference equation 
along the segment 36. Equation (3.14) is used to determine the change 
in S along 86, that is, to find S at 6. Equation (3.12) is used to set up a 
difference equation along 56; the values of /'t’u,; and /'t’’u; ; are taken to 
be those at point 5. Equation (3.13) completes the system, which can then 
be solved by an iterative procedure. When 6 has been found the point 6’ 
is next determined, and so on. 

The process must begin in a surface where conditions are known, such 
as a plane of symmetry. 

In general, it appears that the characteristic coordinate system is 
completely non-orthogonal. Coburn & Dolph made the conjecture 
that the diredtion /‘is normal to ¢' and t’' throughout the flow space provided 
that this is the case on the original initial surface. ‘They did not prove 
the conjecture and in spite of a thorough investigation the present author 
can find no evidence of its validity. However, the obliqueness of the 
coordinate system does not add seriously to the formidable amount of 
numerical work inherent in non-linear problems in three dimensions. 


4. LINEAR PERTURBATIONS OF STEADY SUPERSONIC FLOW PROBLEMS 

At the present stage, when so little is known about three-dimensional 
flow problems, the need for calculating completely non-linear fields of 
How is limited. ‘To begin with, an approximate technique is required 
which, by introducing some simplification in the general method, leads 
fairly quickly to an assessment of three-dimensional effects. ‘This is the 
motive behind the Linearized Method of Characteristics developed by 
Ferri and, in more limited forms, by Ferrari and Guderley. 

We now introduce a Linearized Method of Characteristics from a new 
and more general point of view, starting from the characteristic equations 
set out in §3. Ferri considers the linear, three-dimensional perturbation 
of certain basic non-linear flows involving two independent variables, such 
as plane or axially symmetrical flow past non-slender bodies. Here we shall 
first examine linear perturbations of a general three-dimensional flow and 
then treat perturbations of flow in two variables as special cases. By this 
means it is much easier to ensure that the process of linearization does not 
distort the characteristic coordinate scheme described in § 3. 

We start with a field of flow in which, at every point x’, the velocity 
vector U’, pressure P, density R, and entropy S are assumed to be given. 
We then consider a field of flow which is a linear perturbation of this given 
basic field with velocity vector U'+u', pressure P+ p, density R+p, entropy 
S+s. We assume that the dependent variables in the additional flow field, 
that is, the flow field remaining when the given field is subtracted from the 
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perturbed flow field, are small quantities of the first order, whose squares 
and products can be neglected. We refer both the perturbed and the basic 
fields of flow to the same system of coordinates x', with metric tensor g;,. 





The equations governing the perturbed field are 
(Ui+w)(U,+u,) + — = 0, (4.1) | 
(U*+u*) 
o*(U,+u,),.+ ——, (P+ = 0, 4.2 
5 ( i)k (R +p) +p)(C+c)? , (CF ) x ( ) 
(Ui +w)(S+s); = 0, (4.3) 
P+p = f(R+p, S+s). (4.4) 


In the basic field the same equations are satisfied with all terms in u’, p, p, 
somitted. If we subtract these equations from the corresponding equations 
(4.1) to (4.4) and retain only terms of the first order in wv’, p, p and s we obtain 
the following equations to determine th¢ additional flow field, 


if l ! IT] P = = } 

Uu, + Rpeitu U,;- Re P, = 0, (4.5) 
; l Balt het 7 p 2c -_ = 
gu Da RC 3 l P ,. a RC) u I oS a (i 7 ct Pee = 0, (4.6) 
Uis ,+wS, = 0, (4.7) | 


‘lhe coethcients of the derivatives of uw’, p, and s in (4.5), (4.6) and (4.7) are 
the same as the corresponding coefficients in the equations of the basic 
flow. It follows that the characteristic properties of the additional flow 
are identical with those of the basic flow, so that the characteristic equations 
of the additional flow can be referred to the coordinate system defined in 
§3. When this is done (4.5), (4.6) are replaced by the system 


aay (UN Ct) —a( U4, NCH}, + 
+ aay (ULN-Cl)—a( U4, NCA}, 
” s \/(O? — C?)tip ,- a wP.-(5 + Z)U'P x) 
— aN P,+wnrU,,= 0, (4.9) 
_— (Ut) A" — Ct") — a U1, X"" — Cl) tu, , + 
+ may (UN Cl) —a( V4, CF} Pa, - 


, Caer I : p a 
— REVI CW'P.— olsPa-(G + Z)UPal - 


— FNP, +wNU,, = 0, (4.10) 
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ieee... By ate _ 
+1wU,,— 2.0P,+ 5l'p,=0, (4.11) 


. , > 
(2S nu (4.12) 


Equations (4.9) and (4.10) are the ‘second characteristic conditions’. 
Equation (4.11) corresponds to equation (3.12) for the basic flow, and (4.12) 
is derived from the energy equation. ‘These are to be solved in conjunction 


with (4.7) and (4.8). 


Simplification when basic flow involves two independent variables 

The equations simplify for basic fields which involve only two inde- 
pendent variables, for example, plane flow, flow with axial symmetry and 
conical flow. In all such fields the basic flow is identical in each of a one- 
parameter family of surfaces and it is possible to choose the vectors /' to be 
normal everywhere to the surfaces of flow. With this choice /’ is normal 
to U’, t,t’, X and A’ and a=0. Equations (4.9), (4.10) and (4.11) then 
simplify to 


(Ut, N—Ct')t’u, ,— Clu, ;- Z V/(O?— Cp ; 


~~ RC 
1 ae i, ee 7 
— Beit Pa- (G+ G)U'Pap — faNPatwAU,, = 0, (4.13) 
(AN Creu, — CP, — RVG CNP, 
Bite y) = 
saa RO u” P.-(§ + Z)UP = sat, +4 Nil # = 0, (4.14) 
lic? 





; , p 1 
—— (tu, + tu, ;)+lwU,, — =P;+alp;,=9. (4.15 
naire es > Bet on 
We now develop these equations further in the two special cases of 
basic plane flow and basic axially symmetrical flow. 


Linear perturbations of plane flow 
In the plane flow case we take /' in the direction normal to the flow 
plane with corresponding coordinate z. ‘The characteristic surfaces are 
then cylindrical surfaces and can be written 
% = constant, 8 = constant, 
with generators parallel to the z-direction. Let h,, ; be length parameters 
along the coordinate curves in the «- and f-directions respectively. Then 
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if ~ is the Mach angle of the basic flow, the original coordinate vector x’ is 
Cartesian, with x! and x? in the basic flow plane and x? normal to the flow 
plane. We then find that 
ds? = dxi+dx3+ dx; 
= h® dx? + h3 dB? + 2h, h, cos 2p dudB + dz*. (4.16) 
The coefficients in (4.16) are related to corresponding coefficients 
defined in §3 of Coburn & Dolph’s paper by 
A =k, B = &,, C = {, E = h,h,cos2u. 
Since x' is Cartesian, the directional derivatives of u' become, simply, 
P 1 du. Y 1 cu } Cu; 
he as eg = or tes “4; = =. 
' h; op ¥ h, Cx oz 
If @ is the angle between the basic velocity direction and the x-axis, the 
unit vectors associated with the coordinate system and the velocity vector 
are defined as follows: 
iy icos(#— 2), sin(@— 2), 0}, 
a icos(# + 2), sin(# +2), Of, 
li = §0,0, 1}, 


XN’ = {—sin(@—), +cos(A—p), 0}, 
Ai = {+sin(0+p), —cos(@+ 2), 0}, 
UL’ = {O cos 8, O sin @, O}. 


Making use of these relations, (4.13), (4.14) and (4.15) are transformed, 
after some reduction, to the equations 


Cus 


cot uu op 
h, op Rh; op Oz 





pK+k wt, 449) 


é Ou 
Osin#é —2 ;+ Ocosé 
3 h, op ‘ 


Ou Cus cot 0 0 
O sin @ ea O cos 6 —= — in a 





~C'+KiM=0, (41 














1, O% h, ex R hoa oz 
where 
: cosec2u[, . . . oo 
= = RC | {ut sin(@-+ 4) - u* cos(6+ u ey: - 
f ae 2 ) oP 
+{—u le he Pe a 
. 7 aP 
————— (j= + ia.) (4.19) 
ROsin*u cos \hz 0B hy ea 
\ p 
L = Fa Coty me + cosec 2u{u! sin(@ +) — u? cos(@ + p)} x 
3 OB 
x ¢ —sin(#— jp) ied + cos(@— 2) eg cosec 2u{ — ul sin(@— x) — 
a | f h; op) we ; f 
, dU dU, | 
—u*? cos(@— yz). —sin(@— ~ +-onlb—pi.———>, (4.20 
=p); ~aale— nie tot i 0) 
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C 


M= 





cot u + cosec 2u{u! sin(@ + .) — u® cos(@ + 1)} x 


- | 
"rT 5 
au, aU.) 
21. _. euallts who? 
hap 8" +H) 5 ap, 


p 





>+ cosec 2u{ — u! sin(@— pp) 4 





x < sin(@+ pw) 


aU, 
— —cos(@+ 
hy ex (9+) 


al = al - 1 0 
and TOsec p =m + pa It a2 = 0. (4.22) 
ag h, ep h, 0x * Roz 
The system is completed by the perturbed forms of the energy equations, 
the equation of state and the condition of conservation of entropy in the 
stream direction. ‘These give 


(4.21) 








+ u® cos(@— 2) }< sin(6 + p) : _ rs 
, O% 





_ (pdP (ad 
u,U,+u,U,- Eo { |Z = 0, (4.23) 
oP oP 
P= PaR T8955 oro 
wS ,+Om’s , = 0, (4.25) 


where m’ is a unit vector in the basic stream direction. 


Linear perturbations of flow with axial symmetry 

When the basic flow has axial symmetry it is convenient to refer the basic 
field to cylindrical polars x’ with the x!-axis along the axis of symmetry. 
To calculate the perturbed field we use characteristic surfaces 


% = constant, 8 = constant, 


obtained by rotating the characteristic curves in the basic flow plane about 
the axis of symmetry. Corresponding length parameters h,, h; along these 
curves are introduced. ‘The meridian planes are used as the third set of 
surfaces, so the third coordinate, 


4. 


is the angle between any meridian plane and a fixed meridian plane. In 
terms of the usual notation for cylindrical polars, 
eas = r= ® 

Since the original coordinate space is not Cartesian the directional 
derivatives are not as easy to calculate as in the case of basic plane flow. 

Let y’ be the coordinate vector in the characteristic space with metric 
tensor G.,. Then 

y' = 8B, y = 


ds2 a ae dx” dx” a GCG... dy" dy", 


R 
< 
| 
3 
w 
\l 
c- 


where 
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and a short calculation shows that 
| h= h, h, cos 2u 0 | 
an =| Ay hz cos 2p ‘ 0 
| 0) 0 (ss 
‘The unit vectors in the characteristic coordinate triad are 
t' = {cos(@— ), sin(@—j;.), O}, 


t = |cos(@+ 2), sin(@+ p), 0}, 


G 


The covariant derivative of u, is 


Ou m 
u,; = = 


m). ew: 
where < .. » is the Christoffel symbol. In basic plane flow all components of 
y : 


m 


‘They are 


o>) 


, 
3 
Ke: 23 Lee: 32 1 /x?. 


? 
~ 5) 
> 
3 


» 
by | 


are zero but in the present case three components do not vanish. 


\s a result additional terms are introduced into the expression for the 


directional derivatives along /', t', t’ 


The equations governing the perturbed motion are finally found to 


be as follows: 


Cu Cus Ceu Cun cotmu op : 
O sin 6 — “* O cos @ “wa = : os — +k 
‘ hi,05 * h,cB red } R h;cpB 
; Cu Cus Ceu Ci, cotu cp : 
O sin @— Ocos@ —— —- —— = — ies ae 
F h. ea . h. ez rcd j R hoz 








together with (4.10), (4.24) and (4.25). In (4.26) and (4.27), A, L and MW 


are the expressions defined in (4.19), (4.20) and (4.21). 


Applications 


Ferri has considered a number of practical problems to which his linearized 
method of characteristics may be applied. In the case of perturbed plane 
tlow he takes a finite wing, in any cross-section of which the flow is sub- 


stantially two-dimensional. As an example of perturbed axially symmetrical 


flow he considers bodies of revolution at small angles of yaw and non- 


symmetrical conical bodies. Both classes of problems are important and 
merit futher attention. We shall describe briefly how they may be attacked 
by the present linearized characteristics method. 











Method of characteristics for supersonic rotational flow in three dimensions 423 


In both applications we must first apply the method of characteristics for 
two independent variables to calculate the basic flow. We then know the 
characteristic coordinate network for the additional flow and at each point 
of this we can tabulate the functions K and L and the coefficients of derivatives 
of additional flow variables. The additional flow can then be calculated 
by integrating the three-dimensional equations along the basic network. 
The resulting equations are all linear and the numerical process of solution 
is therefore comparatively simple, although the boundary condition may 
introduce difficulties. Conditions must be satisfied partly on a shock wave 
and partly on the body surface. It is necessary to construct the change in 
shape of shock due to third-dimensional effects as the calculation proceeds. 

Simplifications can be introduced into many problems concerning 
perturbations of plane flow. To the first approximation, plane supersonic 
flow past a curved aerofoil can be treated as that through a simple wave 
behind an attached curved shock with negligible entropy variation. A 
three-dimensional perturbation of this can be treated as a linear non- 
homentropic perturbation of a simple wave. ‘This is governed by equations 
(4.17) to (4.25), where now all the derivatives of the basic variables along 
curves of one characteristic family vanish. 

A similar simplification arises when three-dimensional perturbations of 
two-dimensional wind tunnel flow are examined. Part of the basic flow is 
again a simple wave and it is not difficult to extend the technique of Meyer 
& Holt (1951) to correct wind tunnels for three-dimensional departures 
from specified shape. 


This is the first report describing work carried out in the Dept. of 
Mathematics, Harvard University, under contract N5ori-07634, during the 
period | September 1955 to 31 January 1956. The author wishes to 
thank Garrett Birkhotf for making the visit possible, and for many 


stimulating discussions. 
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SUMMARY 


This paper presents an approximate solution of the problem of 
the onset of convection between plane-parallel plates heated from 
below when the fluid between them absorbs and emits thermal 
radiation. A complete solution to this problem would be 
extremely difficult, and the equation of radiative transfer is there- 
fore developed in two approximate forms, one appropriate to an 
opaque medium, the other to a transparent medium. ‘This 
equation is then combined with the dynamical equations of the 
problem. 

The initial static state is investigated by use of the Milne- 
Eddington approximation, and it is shown that there can be very 
large variations of temperature near to the boundaries. 

The conditions for marginal stability are investigated both 
for motions which are restricted to the temperature boundary 
layer, and for motions which take place in the body of the fluid. 
In the former case it is found that a complete solution is provided 
by the approximate form for a transparent medium, and in the latter 
case a reasonable interpolation has to be made between results for 
the two approximate forms in order to complete the solution. 

‘The effect of radiative transter both on the initial static state 
and on the dynamical equations is such that the fluid is stabilized. 
‘This stabilization could probably be detected in the laboratory 
under certain conditions. In the earth’s atmosphere the critical 
Rayleigh number for large scale motions may be increased by a 
factor 10°, while at the surface of the solar photosphere the factor 
may be as large as 10". 


1. INTRODUCTION 

The gaseous envelope of a star or planet must pass out to space a flux 
of radiative heat. ‘This usually gives rise to motions in the atmosphere, 
so that the atmosphere will not be in radiative equilibrium but will have 
sources and sinks of radiative energy distributed through it. 

Investigators of the earth’s atmosphere always compute the field of 
radiative heating from observed temperature profiles, and attempt to find 
the field of motion consistent with it. ‘This does not go to the heart of the 
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problem, however, for we wish to know what fields of temperature and 
motion are demanded by the radiative boundary conditions, not only 
whether they are consistent with one another. 

The aim of this paper is to explore the possibility of formulating the 
equations of motion and radiative transfer in a way which allows them to 
be solved directly in accordance with the boundary conditions: ‘The 
problem chosen is that of the marginal stability of a fluid between parallel 
plates, heated from below. ‘This is probably the simplest problem of its 
type, and the same problem in the absence of radiative transfer has been 
treated by many writers (see, for example, Pellew & Southwell (1940)). 
The fluid is assumed to be homogeneous in composition and incompressible, 
and the temperature to vary only slightly. ‘The absorption coefficient of 
the fluid is assumed to be the same at all wavelengths, and to be independent 
of the physical state. ‘The upper and lower boundaries to the fluid are 
assumed to be black bodies. 


"THE DYNAMICAL EQUATIONS 
These were treated very fully by Pellew & Southwell (1940), whom we 
shall follow closely. ‘The following symbols will be used: 
v = kinematic viscosity, 
a = coefficient of expansion, 
g = gravitational acceleration, 
z = rectangular coordinates (z vertical and perpendicular to the 








x, VY; 
plates), 
w= velocity component in the z direction, 
ia Sb ee ee, 
ox” oy os” o3° 
6 = temperature, 
ou 

6, = temperature in the initial static state, B = — 
0’ = 6-4, 

H = rate of radiative heating per unit volume of fluid, 
H,, = the same for the initial static state, 
H’ = H-dH,, 

k = molecular thermal diffusivity, 


s = heat content of the fluid per cm? per degree. 
The linearized equation governing the heat transfer by fluid motion in 
a steady state was given by Pellew & Southwell in the form 
—vV4w = agV76’. (1) 
Their derivation of equation (1) made use of the fact that, in the absence of 
radiative transfer, there is constant temperature gradient through the fluid 
in the initial static state. ‘The same equation can, however, be obtained 
with no extra difficulty if the temperature gradient varies with z. 
The condition for a steady distribution of temperature is 


wB = H/s+kV?6, (2) 
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which relates the convective, radiative and diffusive heating; and the initial 
static state is determined by 
0 = H,/s+kV6). (3) 
‘Temperature will be assumed constant over the upper and lower boundaries, 
and therefore H, and 4, are functions of z only. The vector flux of radiative 
energy in the initial static state will be in the z direction, and will also be 
a function of = only. If F, is the zs component of this flux, then 
H, = —dF~./dz, and we may write (3) in the integrated form 
F,—ksB = const. (4) 
Combining (2) and (3), we find 
wB = H'/s+kV6’. (5) 
Pellew & Southwell assume that w and @’ are separable functions of 
x, y and z, and that 


Viw = — Sw, (6) 


where / is the distance between the plates and a is a ‘ characteristic number’ ; 
the same assumptions will be adopted here. 


Further let 
, ey ae 
C¢=({--3 and —= [7 
c=(5 i), =. 


l 
so that V-0 = R (D? — a*)w. (7) 





‘The elimination of 4’ between equations (1) and (5) leads to 
V°w. (3) 


Since f is a function of z only, (8) can be rewritten, with the use of (6) 
and (7), as 








a V?H' = kv (D?-a*bw 
pee=- +2 oe (9) 
We shall seek values of the Rayleigh number 
Baghs 
ag ay eee 0 
oi kv ’ (10) 


where 8 is the mean value of 68 throughout the fluid, at which there is 
marginal stability. In general R will be a function of a?, and:the condition 
for marginal stability is found by minimizing R with respect to this para- 
meter, together with any other variable parameters which are introduced. 

The simplest boundary conditions will be used ; these are that the plates 
are free and conducting surfaces, i.e. 


w= Dw=0=0 at (= +} (11 a) 
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Since temperature is assumed constant over the boundaries, these conditions 
also require that 


D'w =0 at C= +}. (11 b) 


With these boundary conditions the problem without radiative transfer 
(which will be referred to as the conventional Rayleigh problem) can be 
solved exactly, giving the critical Rayleigh number as 


R, = —— = 657. (12) 





3. APPROXIMATE FORMS FOR V7H"' 

‘The function V7H’ can be expressed in the form of a complicated integral 
over the whole fluid. The solution of even a very simple equation of the 
type H, = 0 is a matter of great difficulty (see, for example, Chandrasekhar 
1950), and can only be achieved by iterative methods. In view of the 
present state of mathematical techniques, an exact solution of (9) does not 
seem possible. ‘There are, however, two simple approximate forms for 
YH’, one valid when the fluid is optically thick and the other when it is 
optically thin. When numerical values are inserted, there is remarkably 
little doubt about the way in which these two solutions join up. 

‘The equation of transfer for the problem is (Kourganoff 1953) 


where « is the coefficient of absorption per unit volume, /(s) the intensity of 
radiation in the direction of the vector, s, B the Planck black-body intensity, 
and ds is an infinitesimal displacement in the s direction. ‘lhe radiative 
heating rate is 

| dI(s) 


oS ds 





lw, (14) 


where w is an element of solid angle, and the integral is taken over a solid 
angle of +7. Since the black-body intensity is isotropic, we have from (13) 


and (14) 


wn 
— 


H = —4nxB+k ( I(s) dw. (1 


‘The first term on the right hand side of (15) represents the cooling at a 
point in the fluid due to the emission of thermal radiation at the local tem- 
perature. ‘The second term represents the heating due to absorption of 
radiation emitted by other elements of the fluid and by the boundaries. ‘The 
mean free path of the radiation is «', and the main contribution to this 
second term will come from points spaced at about this distance from the 
point under consideration. From (6) the dimensions of a convection cell 
are seen to be of the order of h/a; so, if x! > é/a, the irradiation at each point 
will originate either far outside the cell which contains the point, or from the 
boundaries. In this case the irradiation will not vary much over distances 
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comparable with a cell dimension, and the variability of H will be chiefly 
due to the variability of B. Hence 
Vit = —4ncViB. (16) 
‘This will be called the ‘transparent’ approximation to equation (15), 


valid when «7h? < a?. 
By Stefan’s law we have that 


7B = of, (17) 


where o is Stefan’s constant. For a small temperature range we may write 
approximately 


oB 4c... 00 oA 
= eee, (18) 
Ox 7 Ox Ie 


and treat O asaconstant. Since 4, and H, are not functions of x and y, we 
have from (18), (16) and (1) that 
ae 47Okv 
Vil ——— V'w, (19) 
ag 
For large x, H can be expanded in a power series in terms of «~?. 
\ formal solution of equation (13) is 


I(s) = ec (’ Ke’ B(a) do, (20) 


where gq remains to be determined from the boundary conditions. ‘lhe 
boundary conditions will have an appreciable effect only at distances less 
than « ' from the boundary, and outside these regions we may neglect the 
contribution from the lower limit of the integral in (20). | Successive partial 
integrations then yield 


4B @B dB 


p—r 1 — 4 peas glia 5 
I(s) B ds ds? ds3 (21) 

and hence using (14), we obtain 

dB d*B ‘BB 
= dwt" | —— dw—«*| —— dw+ ; 22 
H ds d. | ds2 dc K | 7s3 ac re (22) 
Now 

dB cB oB oB e 
ds ity °° ey P8 92? (25) 


where Hy, Me and fig are the three directional cosines of s. B does not depend 
on the direction of s, nor do its partial differentials with respect to x, v, s. 
A typical term in the first integral in (22) is therefore 0B/éx fu, dw, which is 
zero. ‘lhe second integral involves integrals like [u, uw, dw, which are zero 
and integrals like fu; dw, which equal 47. Hence 


2B " 
OF ie = (24) 


|e 3 


| 
| 








ee ee 
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Similar arguments show that the third term on the right hand side of (22) 
is zero, and that the fourth is smaller than the second by a factor of the order 
of a?/(K*h?). ‘Thus, provided «7h? > a’, 

47Ov 


3Kag 


Z 





Vill’ = — Vow. (25) 
We will refer to (25) as the ‘opaque’ approximation. ‘This is essentially 
the equation found by Brunt (1944) in his attempt to form an analogy 
between radiative and conductive heat transfer. The quantity 47Q/(3x) 
has the dimensions of conductivity, and so we introduce the non-dimensional 


ratio 








470 

‘eee 26 
x 3«ks © (26) 

On substituting (18) and (25) into (9), we find that, for «7h? Sa?, 
8 D2 43 3 " 
Rw = = ot ay x)w, (27) 
B a’ : 

and, for «7h? > a?, 
D?— a2)? [ 
Rw! = — ms ) < (D? — a®) — 37h? x > w. (28) 
a 


xe) 


It is of interest to note that so far no use has been made of the radiative 
boundary conditions in deriving (27) and (28), and that these equations are 
equally valid if the boundaries are black bodies or mirrors. ‘The nature of 
the boundaries affects only £, the initial temperature gradient. 


4. ‘THE INITIAL STATIC STATE 
In the initial static state all quantities will be functions of z only, and 
the equation of transfer (13) becomes 
dl 
137. = «[B—T]. (29) 
dz 
Iterative solutions of one-dimensional radiative equilibrium problems 
all show that remarkably accurate results can be obtained by assuming a 
suitable simple form for the angular distribution of intensity. One of the 
simplest assumptions is the Milne~Eddington approximation: 
I(pg,2) = 1,(z), for 0 <p, <1, 
I(pg,2) = I(z), for —1<p,- 


By use of this approximation a simple differential equation can be formu- 
lated for the radiative flux: 


0. Se 


F, = | gl dw. (31) 


On integrating (29) over a solid angle of 47, we find 
dF. 


= = 47«B—27x(I,.+_), ae 
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and, if (29) is first multiplied by ps, 


2rd 
stale | a wr 33 
x 7s [_] KF. (33) 
Elimination of (/., + /_) from (32) and (33) gives 
mic = PI + 3x?F.. (34) 
dz* dz : 


Elimination of B between (4), (18) and (34) gives the differential equation 
of the problem, and we now require the boundary conditions at the upper 
boundary (z=) the molecular diffusion will ensure continuity of 
temperature, and therefore the downward intensity /_ will be equal to 
Bh). Substituting /_ = B(h) in (32), we get 





(Z) ee a (35) 
But, from (31), 
F(h) = x{I,—L], (36) 


and therefore 
dF 
a ee Pees Oy or 
(Ze), ack Ah). (37) 


Similarly, at the lower boundary, 


(=z) = 2x«F(0). (38) 
dz ]o y 


The solution of the problem is now straightforward, and we find 


E = LeoshX+M, | 
B 
2x Ms ? ? . 71 : | 
L=y +. sinh 3A +(3 +3y)!*sinh A+ cosh al :S (39) 


| 

L ; | 

M = —[(3+3x)! sinh $A+cosh $A], 
X 


} 


where 


A? = 3x7h?(1 + y). (40) 


According to (39), 8 8 —-1 if either A or y tend to zero independently. 
If A and y are both greater than unity, there is a boundary layer in which the 
variation of temperature is exponential and which tends to a discontinuity 
as \-- ox. If y >A*, 8/8 becomesa function of A only; figure 1 shows a 
number of profiles for this limiting case. 


5. APPROXIMATE SOLUTIONS 


Pellew & Southwell showed that variational methods can lead to re- 
markably accurate values of the critical Rayleigh number even when the 
precise form of w is not known. The basis of the method is as follows. 
Let the right hand side of either (27) or (28) be written O(w) where O is one 
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of two operators both of which, with the boundary conditions (11), can be 
shown to be Hermitian. Consider the function 


R = =) (41) 
[| w'Q(w') at 
J-4 
where w’ is any function which satisfies the boundary conditions (11). It 
can be shown that R’ is a minimum when w’ is the solution of (27) or (28) 
which leads to the lowest possible Rayleigh number, i.e. the critical Rayleigh 
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Figure 1. 8’B as a function of ¢ for the limiting case y >*. 


number R,. Thus, because of (27) or (28), the minimum of R’ will be equal 
to R.. 

The procedure is to choose a physically reasonable function w’ which 
satisfies the boundary conditions and contains one variable parameter. R’ 
is evaluated and minimized with respect to this variable parameter. 

Two functions have been investigated representing two general types 
of motion which may reasonably be expected to occur. The first, 

w’ = sinnn(¢+4), (42) 
is the correct symmetrical solution of the conventional Rayleigh problem. In 
this case the variable parameter is , which can, however, take only integral 
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values. ‘This function was chosen to represent a disturbance taking place 
through the whole body of the fluid. 

Using the function (42), the integrals in (41) are easily evaluated for 
both the ‘transparent’ and ‘opaque’ approximations. ‘The resulting 
expressions contain both a? and n. ‘The value of a? which makes R’ a 
minimum determines the width of the cells in the least stable mode of motion. 
For the ‘ opaque’ approximation the minimum lies at 





a” = ; (43) 


which is the same result as that obtained in the conventional Rayleigh 
problem. For the ‘transparent approximation’ the minimum of R’ lies 
at a value of a® lying between (43) and 

a? = n*n’. (+4) 
In this approximation, therefore, the cell width is variable, but whenever 
radiative effects are large one finds that (44) is valid. 

When R’ has been minimized with respect to a?, it 1s clear by inspection 
that the smallest value R’ occurs form = 1. ‘The same result is, of course, 
found in the conventional Rayleigh problem. 

Approximate values of R. computed in this way are shown in figure 2 
for y = 10%, 10° and 10°, and for A between 10°! and 10°. ‘The ‘ transparent’ 
and ‘ opaque’ approximations are shown in their respective ranges of validity. 
Discussion of these results will be taken up in the next section. 

The second form chosen for w’ is 

Ww’ (4 + Bl) coth nf + (CC + &) sinh nf, (45) 
where A, B, and C are determined from the boundary conditions (11) and » 
is the variable parameter. If 7» > 1, this expression is very large for values 
of € within about 7! of the boundary values. ‘The reason for adopting this 
function is that the form of 8/8 given by (39) suggests that motions in the 
temperature boundary layer might possibly become unstable before motions 
through the body of the fluid, and therefore that a function representing 
motions concentrated in a boundary layer should be tried. 

It was conjectured that, for small », the critical Rayleigh numbers 
derived from (45) would not differ greatly from those derived from (42), 
since in this case both functions represent motions though the body of the 
uid. Since we are mainly concerned with results differing greatly from 
(42), it was considered to be sufficient to use (45) in its asymptotic form for 
large 7 and to neglect contributions to the integrals in (41) from the neigh- 
bourhood of ¢= 0. ‘This greatly reduces the labour of the computations, 
and it appears to be fully justified by the way in which, in figure 2, the full 
lines and chain dashed and dotted lines converge as A decreases. 

In the half-space ¢ > 0, the asymptotic form of (45) which fulfills the 
boundary conditions (11) is 


pete. o Oe a 2 
w’ = ke” {-}--) dhe (46) 
a 7} 
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Equation (46) and the appropriate operators for the ‘ opaque’ and ‘ trans- 
parent’ approximations can now be substituted in (41), and expressions 
for R’ can be found in terms of a? andy. In doing so, terms of the order of 
x! were neglected in comparison with unity, since the minimum value of x 
used in the computations is 10°. 
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Figure 2. Approximate critical Rayleigh numbers as a function of € for x 107, 10%, 


10°. The full lines, marked ‘ body motions ’, are computed from the function 
(42) for both the ‘ opaque’ and ‘ transparent’ approximations. ‘The dashed 
line is an attempted interpolation over the range where neither approximation 
is valid. The chain dashed and dotted line marked ‘ boundary layer motions ’ 
comes from the use of the function (45) with the ‘ transparent ’ approximation. 


In the case of the ‘opaque’ approximation, it turns out that R’ has a 
minimum with respect to a? but not with respect to 7. ‘This suggests that 
the ‘opaque’ approximation is never valid for these boundary layer motions. 

In the case of the ‘transparent’ approximation, the minimum value of 
R’ was found to lie at A?/n? = 5°11, a?/n? = 0°372, giving as an estimate of the 
critical Rayleigh number 


4r-9 2 
R, = 11-89—[ + = +1]. (47) 


F.M. 2F 
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Computations from (47) are shown in figure 2 for A>10, which implies 
7 > 10(5-11)-"? = 4-4, a value large enough for (46) to be a good approxi- 
mation to (45). 
The ‘transparent’ approximation is valid provided that 
Kh?/a? = A?/{3a2(1+ x)} < 1. 


Since the minimum of R’ lies atA? = 13-7a?, this inequality is always satisfied 
provided yx > 4, so that the transparent approximation provides a complete 
solution for large y. Physically this means that, since the radiation boundary 
layer is always smaller than the radiation mean-free-path (boundary layer 
thickness = hA”! = «-'{3+3y]-!? < «1 for x > 1), motions occurring in it 
can always be treated by the ‘transparent’ approximation. 


DISCUSSION 

Let us first consider the ‘body motion’ curves in figure 2 obtained by 
the use of the function (42). The most remarkable result is that there is 
almost no doubt about the form of the solution in the region where neither 
approximate treatment is valid. ‘This means that a reasonably complete 
and satisfactory solution to the problem has been achieved without entering 
into the usual intricacies associated with radiative transfer problems. There 
is therefore every reason to hope that similar methods used on more realistic 
models will yield results of real importance to the study of planetary and 
stellar atmospheres. 

As far as this particular problem is concerned, the results indicate that 
if either A or x is independently less than unity, then there are no radiative 
effects on the convective motions ; they also indicate that R, = 657, as in the 
conventional Rayleigh problem. If A and y are both greater than unity, 
then the fluid is stabilized. ‘This stabilization has two aspects. First, the 
effect of radiation on the initial static state is to concentrate the temperature 
variations into boundary layers, leaving the body of the fluid in a more stable 
state. Secondly, the radiative transfer tends to damp out any motions 
which may arise by providing a means of heat transfer from hotter to colder 
parts of the fluid in addition to the molecular diffusion. 

The maximum stabilization is achieved in the case of body motions if 
the fluid is opaque, and then the critical Rayleigh number is 657(1+ x). 
This is exactly the result that would be obtained from the conventional 
Rayleigh problem with a diffusivity k(1 + y) in place of k. The same result 
was obtained by Brunt (1944). 

With regard to the ‘boundary layer’ curves in figure 1, the most re- 
markable fact is again that an approximate solution has been obtained by 
comparatively simple methods, and in this case the solution is valid for all 
values of A under consideration. For large values of A, when the ‘ opaque’ 
approximation is valid for the body motions, the critical Rayleigh numbers 
are much higher in the case of the boundary layer motions, and so these 
motions should not occur. For smaller values of A, the two sets of curves 
converge, and in general the boundary layer motions are more stable except 
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near A = 10? for y = 10%. It is very doubtful whether much physical signi- 
ticance should be attached to the apparent change of regime that the mathe- 
matical analysis indicates. It is much more likely that, where the boundary 
layer and body motion curves lie close to each other, the motion is very much 
more complicated than equations (42) and (45) would suggest, and that they 
both give approximate critical Rayleigh numbers which are considerably 
too high. 

The results presented here may have some bearing on a suitable labora- 
tory experiment. Large radiative effects are more likely if a gas rather 
than a liquid is used as a fluid; but, unfortunately, radiative transfer in a gas 
at normal temperatures will be by vibration—rotation bands of great com- 
plexity for which a mean absorption coefficient has little meaning. Never- 
theless, a crude estimate can be made for water vapour at S.T.P. from 
data given by Cowling (1950); this is « = 2x 10-*cm™! (the uncertainty 
in this figure may be as high as an order of magnitude). With 
7O = 6x 10% erg cm “sec”! deg"! (corresponding to 273 K) and ks (thermal 
conductivity) = 2 x 10% ergcm™! sec"! deg ', there results y = 2 x 10? and, if 
h = 10cm, 4 = 3. Inspection of figure 2 suggests that the critical Rayleigh 
number might be twice as great as without radiative transfer, which is a 
ditference which should be detected easily in the laboratory. 

Planetary and stellar atmospheres have been mentioned, and it is of 
interest to record the magnitudes involved. In the earth’s atmosphere 
x might be of the order of 10°, and motions on the scale of 1km can be 
treated with the ‘opaque’ approximation. Just beneath the surface of the 
solar photosphere (van der Hulst 1953, Minnaert 1953), y is of the order 
of 10!*, and motions on the same scale as a solar granule must be treated 
with the ‘opaque’ approximation. 


The author is indebted to Dr W. V. R. Malkus for many interesting 
discussions on convection problems, and to Mr J. 5. A. Green for checking 
this paper carefully. 
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SUMMARY 
Propagation of spherical shock waves through self-gravitating 
polytropic gas spheres such as stars, caused by an instantaneous 
central explosion of finite energy EF, is discussed theoretically. 
The problem is characterized by two lengths Ry, L, where 


7 \1/3 3C2 \v2 
Kk x (=) : ZL £8 (; 5 .) ; 
4p 2mpyG 


Po» Po and Cy are the values of pressure, density and velocity of 
sound at the centre of the equilibrium pre-explosion state, and G is 
the constant of gravitation. , and L are scales connected with the 
power of the explosion and the dimensions of the star respectively, 
and their ratio A = R,/L has a fundamental significance. A 
solution especially suitable in the case of A = O(1) is developed 
in the form of power series in R/R, (R is the distance between the 
shock front and the centre) by a method similar to that used in 
previous papers by the present author (1953, 1954). An approxi- 
mation to this solution is carried out up to the term in R*. In 
particular, the velocity of the shock wave U is found to be 


U R \-32 R\2 R \3 
— « 60a <1+0°41.47( —} +0°57[/ —) +...5 
‘ol (z) (x) (Fe) J 


for the case of y = 1°4, where y is the ratio of specific heats. 





1. INTRODUCTION 

‘The purpose of the present paper is to investigate the propagation of 
spherical shock waves through gravitating gas spheres, such as stars, 
caused by an instantaneous central explosion of finite energy. ‘The problem 
has previously been treated by Kopal (1954) and Sedov (1954) for the case 
of constant shock front strength in a special distribution of density in the 
pre-explosion state. The case of variable shock strength is discussed by 
Lidov (1955) and Rogers (1956). We consider the problem more generally, 
with variable shock strength and the pre-explosion state simply given by the 
equilibrium equation of self-gravitating polytropic gas spheres. In the 
case of uniform pre-explosion states, the author (1953, 1954) discussed the 
problem by a method of constructing the solution in power series of U~?, 
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where U is the propagation velocity of the shock front. A similar method 
is also applicable to the present problem, if we use R (the distance between 
the shock front and the centre) in place of U-?. We cannot now use U~? as 
an independent variable as was done in the previous papers, because U is 
no longer necessarily a monotonically decreasing function of increasing R, 
but may even increase as a result of decreasing pressure, density, etc. in the 
equilibrium state, and would then cease to be single-valued. An approxi- 
mation to this type of solution is carried out up to the term in R®. The 
terms in R° and R* can be found from the results of the previous papers, 
since they contain no effect of non-uniformity in the pre-explosion state. 
Further, the term in R vanishes automatically, and the effect of non-uniform 
pre-explosion distribution comes only from the term in R? in this approxi- 
mation. Since U~? = O(R*), this just corresponds to the second approxi- 
mation in the previous paper (1954) including terms up to U~?. 
‘The present problem is characterized by two lengths Ry, L, where 


E \i3 
R, — (=) ’ (1) 


3Cz2 \12 
a 2 
b= (sa) 2) 


and E is the energy of explosion, py, py and Cy are values of pressure, density 
and velocity of sound at the centre in the equilibrium pre-explosion state, 
and G is the constant of gravitation. R, represents a scale connected with 
the effective range of the power of the explosion and appears also in the 
previous papers (1953, 1954). ZL represents, on the other hand, a length 
connected with the dimensions of the star. ‘The ratio A of these wwo lengths 
detined as 


7 


13 
A= R, L, or A? = (4) Po GC, "Po a (3) 


has a fundamental significance. The present theory is specially suitable 
for the case of 4 = O(1), that is, the case in which the effective range of the 
explosion is of the same order of magnitude as the scale of the star. In the 
case of a weak explosion (A < 1), the theory is still valid but reduces to that 
for the case of the uniform pre-explosion state; it then seems to be more 
suitable to consider the problem by Whitham’s method (1953). In the 
case of A > 1, we havea very strong explosion, far bigger than the star itself, 
and it is clear that this phenomenon requires completely different formu- 
lation. 

We formulate in §2 the fundamental system of equations given by 
equations of motion, conservation of the total explosion energy and boundary 
conditions of the shock front. ‘The solution in power series of R is developed 
in §3, and the term in R?, which is the only new one required for the fourth 
approximation, is found by numerical integration in §4. In §5 and §6, 
miscellaneous results obtained from this fourth approximation to the 
solution are discussed, 
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2, FUNDAMENTAL EQUATIONS 
The equations governing the spherically symmetrical flow of a poly- 
tropic gas of adiabatic index y under the influence of its own gravitation 


are (Kopal & Lin 1951) 


Du lc Gin 
= 4 ne (+) 
Dt pcr r= 
Dp (cu 2u\ A 
) F\ - T ar (>) 
Dt er r 
“mM ; 
= 4rro, (60) 
. 
Dy pp = 
5 ile (7) 
Dt 
here m(r,t), u(r, 4), pv, £) and p(r, t) denote respectively the mass inside 
sp! ere ol d . | velocity, Cire pre ssure and the density ata distance } 
trom the Orig at tink t,and the expr sion D dD: denotes D Dt C ct uc 5) 
By means « kation (5), equation (7) transforms into 
Dp ul es 
faa (NS) 
Dt "\ oF } ) 
bor the equ liby um state ' ' / ) nf O. we put 
p p (R), p p (RR), Hl wm (RR), } aK: )) 
thus equations (4), (6) and (7) are reduced to 
l dp Gm dm ; ; 
r : (), ; 4- e205 ’ b'p Po p : 14)) 
p dR R: dR ; 


while equations (5) and (8) give merely the relations p = p’ and p Pp 
this case. 

For the radius of the shock front it will be convenient to use the same 
symbol R as above. ‘This R is a function of ¢ and connected with the pro- 
pagation velocity at the front U by the relation 


dR/dt = U. (11) 


At the shock front (r = R), we have the following Rankine-Hugoniot 


conditions: 


u 2 C\2 p 2y (U\2 y—-1 
U- y+i , (=) i p’ ate ) eae 


p y+1 : Cy) 
aes 4 1 + ome B) 2 m= a, 
p dio — 


where C is the sound velocity in the equilibrium state and is given by 


ayfy’\12 
c- (4). (13) 
p 
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Since we consider the case of an instantaneous explosion of energy E, 
the equation of energy may be written 


-R Ec ’ ~ , 
| (5u+ alight « me 4a? dr—| (S I )ptrR® dR a K 





ott y—lp r —1 p’ R 
(14) 
In place of r, t, we now introduce new independent variables x, 3 detined 
by 
r R S 
— = X, — = 2, 15 
R R, ia 


and express the quantities u, p, p, m; p’, p’, m’ as follows: 


' U\? ) 
7s f(x, 2) -_ P(Z) 2(x, 3), f = Dh(x, 2), - = Mi(x,z), (16) 


U Pi Po 0 
p’ p’ m ; ye 
— = P(z), — = Jz), — = M(z), m, = =p, R3, (17) 
Po ) Po Mo ) in 


where f, g, h, i; P, D, M are non-dimensional new variables. Using (15), 
(16) and (17), we have 


a_1a D Wig et ss 
dr = R ax’ Dt — R\S~ 9 5y +7 5a" a 
We now substitute (15), (16), (17) and (18) into the pane equa- 
tions (4), (5), (6), (8), (10), (12) and (14). ‘Then (4), (5), (6) and (8) become 
_ Sf Ff a2dU, . le , a ony 
30 4*32+ Dae! yas Wee | MY) 
. .loh z0h sdD_ of a ‘ 
UD nay + hae t D dz ~~ ax 7 x? it 
ai ya D 
7 32 we h, (21) 
. og z= og z dD 22 dl of 2f 
G-9, 2a’ ba’ ta” -(2+2), (22) 


while equations (10) become 





1 dP . 2yA? dM 
D Py + pa M= U, > = 32? D, PD’ = I. (23) 
where we have used (1), (2) and (3). The boundary conditions (12) become 
a 2 -1/C\? 
( —! Pe. Oe a 
filns)= —{1-(F) bale) = 2-225). 
- “ 1/ 2 /Cc a , " 


while the energy equation (14) gives 


jz g 
: ci >| G hf?+ — i dx —2yA®MD2* | hix dx — 
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sais 
If we write = | (Sar f) dx = J, 





os 2 MD 
| (—, —2yA? —)e 2dz = I, 
“0 ; iam 
— U? | 
our equation is reduced to ae Dz] —2yA®?MD2K — I = 1, | 
ive + = (25) 
or G2 = PID +1+2y4*MD2K)". | 
J 


> 
~ 


‘THE SOLUTION IN POWER SERIES OF 
Following the method used in the earlier paper (1953), we construct the 
At first we find the solution for the equi- 


Eliminating M and P from (23), we have 


solution in power series of 2. 
librium state expressed by (23). 


, 4D 


“ (26) 


)+o4 2D = 0), 


from which we get 
13 — 
D = 1—A?z? 4+ —— 


(27) 


13—5y tok 2 
10 ANs 


and then, from (23) and (13), P, M and C are given by 


P = 1-—yA?z? + 2yAt24+ ..., (28) 
M = 23 1—34?32+ sas (29) 
C2 ae , ae 2 

C2 = 1 —( = 1) A?2? + Fy - 1)A - (30) 


In figure 1, the expression (27) for the case y = 1°4 is compared with the 
exact solution of (26) which was given by Eddington (1930). 
We now assume 








(x, 2 


f( 
9( 


x,2) = 


where f,(x), g,(x), h(x), i,(x) 
determined. 


We insert expressions (27), (28 


and then obtain 


J = J((1+ 24, + 272, + 23x53 +...), 


) = fot fi tart zit. 


Sot 281 
= hyh+zh,+ 


+ 379,+2 3g, + (31) 
shy + hy + sei 


X, 3) = lg +3, +3%g+2igt.., | 


(v = 0, 1, 2, ...) are unknown functions to be 
), (29), (30) and (31) in equation (25), 


(32) 


ie 
= 
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where 





SN 
ll 


(z 5 hy FZ Bo _\2 dy, 
/0 , 1 


1 
%Jq = | (3 fo+rfolo fir ft |e dx, 
~0 se 





teJy = k 5 he fi+rvfhoMo fe hate cial dx, ¢ (33) 


tg) q = " {asf + foo f+ hal (h, fit 2fy Az) 


Jo 








+Y(fofoiithoh fo) + = i +x? dx, 








J 
1.O 
2 
O. ST 
| 3 
{ { 
O O.5 1.0 Az 
Figure 1. Comparison of approximate density distributions for y = 1:4 with exact 
one. Curve 1, 1—(Az)?; curve 2, 1—(Az)?+0-6(Az)*; curve 3, exact 
solution. 
and K = K,(1+2B,+2°By+2°B5+...), (34) 
where : Ee, ) 
Ky = | Agigx dx, 
“0 
a | 
By Ky = | (hyigthoi,)x dx, > (35) 
~0 | 


rl 
B. Ky = | (hy iy + Ro in +h, i;)x dx, 
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and MD = 2*°(1 —;A?2? +...), (36) 
Yi) 2) os. 3) ied 37 


and hence 





z dU ae 3 1 1 
ag me 3 de we —(e,— A dati ie — 5 + 
U dz 2 =r _ —— i. oad 
b4(4; — 3ay) »39 
= —-(1+6, Og 3° +05 2 ) 
where 
6, = 4a, 5. = 3(2.—A?—4a?), 8, = a3- = me 1 (a? — 3a9). 
ay = 


(39) 


Equation (38) gives the relation between the propagation velocity U and the 
position R = R, = of the shock front, if we know the values of constants 
J) 1, Xe, «. for a given A. 

Now, substituting (27), (28), (29), (30), (31), (38) and (39) into (19), (20), 
(21) and (22), and comparing the coefficients of the same powers of z on both 
sides of (19), (20), (21) and (22), we get the following system of equations: 


for the term independant of 


, £0 
(fo—*)\fo— sho he’ 
hi ee 
(Jo—) h, = a xo ‘ (40) 
go! Pay 
(fo-—*)= = 3- = fo-vfo» 
ly = 3x7hy ; J 


for the first power of = 


f I , ¥ 9 a } 
(fo fi + 81 = B-SMhit Sahat Body | 
ae = Oe | 
Uj) t= -(G + 5)h- (41) 
ow _(8 4 27\,_& 
(fo 0)(&) yf = (@ Yh + 36), 
1, = 3x%h,; J 
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for the second power of = 


; " oe nc —— — 
(fo—x)f2+ 82 = —(3+fo)fe- 2a +5 fo 24 


yho 





a h,[., Ay 
+fi( —J, +50,) ss sae 50 Z) 


. h,\’ /(k,\'. . (5 
eR) (orl) + 02 


\f82\’ r &o  2y 782 942,25 | 
(fo- “‘_: + yf, SS - fz 4— + 2A? + 36,+ 


Oo 
5 0 


o,/o ee? g,\= 
1 1 é 3 $1 
(fo— x)= A = i) f+(%) ’ 
£0\8o S50 £0 


for the third power of z 


] : ree g. 
(To v)f., + vite 23 _ is t fos t me h., r 503 fo t 


g—X + f = fig~ 2 + 
hifh\' (hh? hy\(hy\’ (7) h,(h,\' 
“(lo 9) Rae) ~ (ig Ba) FAA) — BARE) f= | ™ 
-s{ts) — (4a) 
a hy \he he} 


” 3 pees ad 
&i\ | ¢(/82\ — £:f8\’ 
oa Nie od eee od 
8o/ | 80 S0\8o/ J 


oO fa-\* ge a. \ 
+(h- 1), (2) -(& - &) 
Jo §80\S8o 8 8&0 
sf) -9 (8) 
£0 80 \Z0 §o 


Ah) 
= J 
where the primes (’) denote differentations with respect to x. 
From the boundary conditions (24), we have in a similar manner 
2 2, +1 
a 2y y ; 
’ 1 —_ ’ h 1 Se 1 => Pp 
yp = 5, ay", t=, (4m 


~ 
II 
Ww 
= 
nw 
= 
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f{i1) = 0, g(1) = 0, h(i) = 0, (1) = 0, (45) 
Jil) = 0, g(1) = 0, h1) = 0, t(1) = 0, (46) 
2 P ee _  Ay+)) - 
f,(1) sd >+1” (1) a ae ea h,(1) aie Sie Ga1pP” (1) = 0 
(47) 


‘The first step in the solution of the problem is to solve the system of 
differential equations (40) about fo(%), go(*), A(x), io(x) with the boundary 
conditions (44), and to insert this solution in the first equation of (33). 
Then we have the first approximation to the solution in the form 

- U\2 : 
u= Uf(x), p= r( 7) £o(*), Pp = Poho(x), mm = myig(x), 


Co ‘ ‘ 
( U ) Si ad (48) 


Since i, enters only in the fourth equation in (40), we can find fo, go, Ao 
separately from the first three equations. It is clear that this procedure is 
precisely the same as that of the first approximation in the case of a 
uniform pre-explosion state, which was first developed by Taylor (1950). 
The solutions f(x), go(x), Ap(x), 7») were, in the previous paper (1953), 
written as f(x), g(x), A(x), Jo. The value of J, for y = 1-4 is 0-596, 
and some values of f), gy and A, are shown in table 1. 

‘The second step is to find the functions f,(x), g(), A,(x), 7,(~) and a 
constant 6,. For this purpose we use the system of differential equations 
(41), with the boundary conditions (45) and the integral condition (the 
second equation of (33)) with the relation 6, = 4a, from (39). The equations 
(41) and the integral condition are, however, linear and homogeneous, about 
fis 21» Ay, ty, 5;, and so we have simply 


f(x) = 90, g(x) = 0, A(x) = 90, (x) 0, 8,=2%,=0. (49) 


Equations (49) lead to a simplification of the equation for the third step, 
and we get from (42) and (46) 





, I jen l ; £o 3 
(fo—x)f yh,?? = -(5 “fi \h + ae ) 540 Oo, 
h,\' — 2 h 
( — ef el =i ae 24242 
(fe »( a) ri ( 2+ <) fe i +2 
og , ao’ ? o a 
(fr-3)(22) + vfo = -(& “3 2 \f,— 2% 2 Az4 35,, (50) 
So Lo £o 
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and from (33) and (39) 
<a 
tJ = |, (vt fo+vFo%ofe+ 





—8)s* dx, (51) 


82 = 3(a.—A?). (52) 


By means of (50) and (52), we can find the functions f,(x), g(x), Ao(x), to(x) 
and the constants a, 6,. The solution is given below in § 4. 
Equations for the fourth step are also simplified by means of (49), and 
we get from (43) and = 
es a 3 
(fo f+ 58h = - (5 +f fat Bohs + 585h 


P is | = 





j hs = h,, 2 ahs 
(fo- a(t) a = -(j Ls =) fs o % | 
g) 2 a ali ae | 
(fo-3)() ai vf’ = - -(2 a 2) f,—38 + 303, (53) 
50 : 50 
1, = 3x*hg, 
2 , yyt+1). | 
Y ay ) | 
f,(1) = “7 g3(1) = er h3(1) = - (y—1)? ail 
il) = 0, 


and from (33) and (39) 
1 /y : 
tgJy = | (3 hs fo t yhy fo yi 

7 \@ 





1 
5, gs)" dx, 





. 1 
Os = %&— om i 
Now we put 

fs = J o(x —fo)ds; 83 = Jo Bos; hs = Joho X3) d3 = JT Az. (55) 


Then we get 





oe ae re ee 
— (x ~— e (2 wares 5 x) | 











f 
— #3) + 5 5 os 
(x —fo)( — $3 + x3) = 35+ X3)s 
(x—fo)(—vb3+ 453) = 3i(y—1)bs + os —Ag}, + (56) 
i, = 3S q x*hy Xz, 
z y-1 ns 
$3(1) = rrveg 1 #1) = — “2y ’ xa(1) = — eas 1,(1) = 0, 
(" Lo - Ls re 4g 
shy h o(* —fo)bs + ine 5 fio X3 px? dx = AgJg+ 2 f 
hi P ) ae J 


where we have used equations (40) for fo, go, Ao. 
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By comparing this cena of equations (56) for $3, 3, x3, 23, Ag with that 
for 4, #, x, A, expressed by (2), (3), (4 ) (5) and (6) in the previous paper (1954), 
we can find that in the case o a = 2 (spherical wave) they are precisely the 
same except for the equation for z,. ‘Then we get 

$3(x) —- P(x), s(x) = w(x), xa(*) = x(%); As A). (57) 
The value of A, for y = 1-4 is — 1-918, and some values of $3, ob, y, are also 
given in table 1. 


\ fo(x) Lo(X) ho(x) h(x) Us3(x) x3(x) 
1-00 0-833 1-167 6-000 5-000 —(Q-143 5-000 
0-95 0-751 0-760 2-464 3-941 0-779 0-879 
0-90 0-685 0:595 1-234 3-603 0°475 2°567 
0-80 0-584 0-471 0-392 3°45 0-18 3°32 
0-60 0-429 0-427 0-041 0-51 
0-00 0-000 0-426 0-000 3°5 0°53 335 


Table 1 


We have now arrived at the fourth approximation to the solution, 
including the terms up to 3°. Since U-? ~ 33(1+...) by (38), this fourth 
approximation corresponds to the second approximation of the previous 
paper (1954); because, in the former case, the solution was expanded in 
powers of y = (C/U), where C denoted the velocity of sound in the uni- 
form atmosphere. ‘The fact that we can express the coefficient of 3° in the 
present solution in terms of the coefficient of y in the previous result is due 
to this reason. But it must be remarked that the errors involved in the two 
approximations are not the same, since the error in the present case is 
O(2*), while in the previous case it is O(y?), i.e. O(2*). ‘Thus, to have the 
same error, we must proceed to the sixth approximation. ‘The non-uniform 
pre-explosion state, caused by gr werwe and expressed by the constant 4, 
does not affect the terms in z° and 3%, but appears separately in the term of 
3? in the fourth approximation. In the case of A = 0, equations (50) and 
(52) become homogeneous, and the term in 2 vanishes as in the case of the 
term in z; the solution then reduces to that in the case of a uniform pre- 
explosion state. 


4. SOLUTIONS FOR fy, £0, ha, 5s 
We now solve equations (50) and (52) to the third approximation. ‘The 
procedure is similar to that of the second approximation (the term in y) 
developed in the ee paper (1954). We put 
fy = A?(x—fo)do, 82 = A’go te, hy = AN x2 at, = A’o,, (58) 


so that (50) and (52) are reduced to 





x — f,)bo + So ws 
Joe yhy(x =n)" P 


fo 





= (}-2f))d.+4 a — a) + og—1), (59 
(1 2fiiba + IE (tam va) + Se (42-1 (59) 
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(x—fo)(—$2+ x2) = 3¢2+2x2-2, (60) 

(x — fo) — v2 + 42) = 3(y — 12 + 22 — 209, (61) 

iy = 3Ahy x*(x2— 2), (62) 

_ = (1) = xa(1) = (1) = 9, (63) 

oe ee {4 fy ho( x — fp)d: re hot Shih xe -3? dx, | a 
= cae — 1), 


where we have used (40). Since 7, appears only in (62), we can exclude it. 
The problem is then to solve the system of differential equations (59), (60) 
and (61) for ¢., 2, x2 with the boundary conditions (63), and to put this 
solution, which includes o,, into equation (64) and determine o,. It is 
remarkable that A has no role in this process. 

Just as in the previous papers (1953, 1954), we have an intermediate 
integral from equations (60) and (61). Let us perform the operation 
(5y — 3) x (60) —5 = (61); we have 

3h — Shp + (Sy — 3) xo = 2{3h2 — Sie + (Sy — 3) x2 + 502 — Sy + 3}. 
Integrating this and inserting the condition (63) to determine the integration 
constant, we obtain the relation 


. “ 2dx - 
3h. — 5g + (Sy — 3) x9 + 50n—5y +3 = (50,—5y + 3) exp( | = ) (65) 


\- 1 x —fo 





Eliminating 43 from (59) and (61), we have 
bo = ads + bibs + Cx2 + ay + d, (66) 
where 


3 y—-—1)\, . . 
| a 3¥ £25) al 
a (2f a+. Nf d x—fo f 
—s oe (fe, +2) 
b = > ot So se 6 ’ — Pges jee we ’ 
(7 2 x—fo ay : x-fyo y¥ J 
To. ( F : XJ f (x—fo)° 
— (4+ 2 3b a aa ho ian —— " 














50 \ 
Now we split up ¢5, #2, x2 into two parts 
$2 = $21 +2 do2, hy = Poy + Fo foo, X2 = X21 T F2 X22 (67) 
and substitute these in (66), (61), (65), (63), (64). ‘The equations split up 
into two independent systems for $51, 21, X21 ANd dy0, Yo2, Xoo respectively, 








viz. 

ho, = Aba + bbe + CX +4, 

: Fr 1 
bi, = 2] ¥- 7 BO- Mbt 2a || 
" <a (68) 
me “© 2dx 

Xa = Sy a3 en — Sia) + 1 - exp | x =) 

with dull) = %(1) = xa(1) = 9, J 
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and 


with 


Moreover 


where 


I: or 4 
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hoo = Ado» is a bibos Bs a CX29 7 é, | U | 
. po = : } tal 
I -_ . ™ 
$22 S| vie —= 13(v — 1) b22 + 222-2} |, 
Y —fo { (69) 
ena (* 22) 52+ Sh 
Xee = £4) XP : — I~ IPo9 + Yo2 0s 
sy—-3 | J1*¥—fo/ J 
doa(1) = fo0(1) = xXo0(1) = 0. J 
doJ J, + dyJ 5, Or Gy = J, (J) —J.); (70) 
-] ) 
3 4 P41) - 2 - | 
J, = | 4 fo M(x —fo)ber + ci 1 fo, + 5fo hy X21 ‘oll dx, 
/0 - J 
Je =| df hol*—fobert 2 Y £2 hay X09 bx? dv 
J, = | 4 ¥fo h(x —fo)be2+ y-1°2" afi 49 X22 (XX. 
i oad 7 
- 1-4,(68) and (69) were integrated numerically from x = 1 with the 
initial values $y, = 21 = Xo = $o2 = Poe = X22 = 0. Fromx = 1tox = 0-96, 
the Runge-Kutta method was used, taking the steps of the numerical 
integration as 0-02. In the remaining part of the range, the Levy method 
was applied, again with a step of 002. Inserting these values of ¢.,, a), 
X21; doe, Woo, Xoo into (70), we get the value 
co, = 0-182, with J, = 0-596, A, = —1-918. (71) 
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Figure 2. The solutions ¢o, te, x2 for y = 14. 
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Using this value of o,, dy, %, yo were determined finally by (67); they are 


tabulated in table 2, and are shown in figure 2. 
4 


X fo(X) by(x) X2(%) 
1-00 0-000 0-000 0-000 
0-98 0-037 0-081 0-232 
0-96 0-061 0-124 0-386 
0-94 0-076 0-146 0-495 
0-92 0-086 0-157 0-572 
0-90 0-091 0-159 0-629 
0-88 0-094 0-158 0-671 
0-86 0-097 0-155 0-704 
0-84 0-097 0-151 0-730 
()-82 0-098 0-147 0-752 
0-80 0-098 0-143 0-764 
0-78 0-10 0-14 0-78 
0-76 0-10 0-14 0-80 
0-74 0-10 0-14 O-81 
0-72 0-13 0-81 
0-70 0-13 0-82 
0-68 0-13 O-$2 
0-66 0-82 


‘Table 2 


5. VELOCITY—DISTANCE CURVES AND TIME—DISTANCE CURVES 


Inserting (49), (55), (57) and (58) in equation (38), we get 


By use of the values /,, a, As given in (71), and with R/R, in place of z, this 
equation becomes, tor y = 1-4, 


Co . a | R iy RD A2 R ‘ " R : 
( a) = ()-: 6( R) | | —) S2A ta — | (ze) coe 7 
( R\-32 2( =) s7( R\3 = 
a Bae + () 2( —- } +0-57[ —} ...>. 72 
or A 1-30( =) 140-414 Re) 0. Re) (72) 


0 0 \ \ 
Velocity—distance curves (U-R curves) given by (72) are shown in figure 3 
for A* = 0, 1, 2, 5 and 10. 

We can see in figure 3 that 4 has the etfect of preventing U from de- 
creasing, and this effect increases with larger 4 as expected. It is interesting 


to consider UC, the strength of the shock wave, which is obtained by 
multiplication of (72) and (30). In the case of y = 1-4, we have 


C\2 Cy\?/C \2 R \3 R\2 R\3 
: \ sd “0 ()-596{ — ; = Rae ( =, * ( ‘) 
(ic) (7’) (&) ‘ (z,) i ia a 


[ i R \ 32 r R\ col R\3 on 
or a = -30( 3) ? + 0-614 (x) +057 5 J ve (73) 
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from which we can see that 4 is more effective in maintaining the strength 
of the shock wave. 

Using the expressions (73) and (12), we can deduce relations between 
the radius of the shock front R and quantities such as pressure and density 
at the shock front. ‘lo obtain the most interesting of these relations, we 
consider the temperature 7' at the shock front; we have 


L pp _ wy a *T ( 2 y-1 (5 z 

r p’ p (y +1) 2) Y- l 2y ) l ) st lke 
where 7” denotes the value in the equilibrium state. Since 7”/7,, = (C/C,)°, 
where 7,, denotes the value of 7” at the centre, we have finally 


z 0:326( z ) : l -us2ae(e) 4042) woos (74) 





20: 
—§ 

= 

fe) 
0-2 0.4 6 O.8 R 
g2—— 
Ro 

Figure 3. Velocity—distance curves for A? 0, 1, 2, 5 and 10 
In this case, 4 does not have so much effect. ‘The distance—time relation 


can be obtained from an integration of (72) with use of the relation (11). 
This gives 


and we get 


Cy R\°? ‘R\ R\3 
0-308 . 0-23.42( 0-27( Laas 75 
R, s( R ) R ) R ) ae 
42 8 
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Figure 4+. Time-—distance curves for 4? 0,1, 2, 5 and 10 
If we combine (75) and (74), we can get a relation between the temperature 
T and the time f. 
6. DISTRIBUTION OF VELOCITY, DENSITY AND TEMPERATURE 
The velocity u, pressure p, density p and temperature 4 behind the 
shoek front are given by (16), (31), (27), (28), (72), (73), (49), (55) and (58), 
that is, - ( 
C= Eur fo)( A? 2*b2 + 0°5962"43...)§, 
; (76) 
... 1-30¢ 2(] 0-41 4222 + 0-573 ) 
( ia} 
p U 
P{ —]} 2,{1 + A227, + 0596275, +...}, 
n p ( “ 
). 7" ‘ee 
P 1 —1-4A%e?+..., (77) 
U\2 ae ea 
e) = 1-682-3/1 + 1-22-4222 + 1-143? +...}, 
— = Dhy\1+ A®s?y, +0-596z%y_+...4, ™ 
Do X2 peers (78) 
) D = 1—A?2?,,,, 
; 
me se a, (79) 
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\ special interest attaches to the variation of p and p with time, as shown in 
figures 5 and 6 for the case of A* = 2. In figure 5 distributions of p/p, 
when = = R/R, = 0'1, 0°3, 0°5, which correspond to Cot/Ry = 0°001, 0°014, 
0046 respectively, are shown with the initial pre-explosion distribution 
D = p' p» = 1-22". We can see that p/p, increases to several times the 
initial value D at the shock front, and decreases to zero toward the centre. 
‘This decrease is so rapid that there is always a region of almost vanishing 
density near the centre, and this region is also expanding with time. ‘The 
same kind of distributions of 2*p/p, for the same values of z and the initial 
distribution P = p’/p, = 1—2°82? are shown in figure 6. In this case, the 
increase of p p, from P at the shock front is very large, especially for small z, 
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Figure 5. Distributions of density for 42 — 2 at various stages. Curve 1, initial pre- 
explosion distribution D; curve 2, distribution of p/pp) at z R/R, os I 
curve 3, distribution of p py at 2 R/R, 0-3; curve 4, distribution of p pp 
at 2 RR 0-5 
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Figure 6. Distributions of pressure for 4? ~ 2 at various stages. Curve 1, initial pre- 
explosion distribution P; curve 2, distribution of z°p ‘po at 2 R/R,=G6'1; 
curve 3, distribution of 2°) py) at 2 R'R, — 0.3; curve 4, distribution of 


s*p/po at 2 R/R, 0-5 
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and is about two thousand for : = 0-1. The value of *p/p, is shown in 
the figure, though P is depicted without change. Corresponding to the 
region of vanishingly small density, we now have the constant pressure 
region of about the same size. ‘This value of constant pressure changes 
rapidly with z, being roughly proportional to z~8 or (U/C). The distri- 
bution of temperature, which can be obtained from (79), using the above 
values of Pp py and p/p», then has a region of almost infinite temperature 


corresponding to the above regions. 
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REVIEWS 
Fluidization, edited by D. F. OrHMeR. New York: Reinhold Publishing 

Corporation, 1956. 231 pp. 56s. or $7.00. 

I luidization is a modern chemical engineering technique which raises a 
number of new problems in fluid mechanics. ‘The word describes a process 
whereby a bed of small solid particles is given the properties of a fluid by 
blowing gas upwards through the interstices, and may best be introduced by 
describing a typical laboratory experiment. Suppose a vertical pipe con- 
tains sand, which is supported to a depth of a few inches by a wire gauze, 
and that air is blown up through the bed. At small rates of flow, the air 
stream suffers a pressure drop depending on the flow and on the size of the 
sand particles, but does not affect the arrangement of the bed. Ata sufficient 
How rate, however, the overall pressure drop is enough to support the 
weight of the bed, which must then expand to accommodate any further 
increase in flow. ‘lhe expansion may be as much as 30°,, and the fully 
expanded or ‘ fluidized’ bed has many properties in common with a liquid; 
its upper ‘surface’ remains horizontal when the containing vessel is tilted, 
and it offers little resistance to the movement of objects floated on the 
surface. If air is blown at a still greater rate through the fluidized bed, it 
ceases to escape uniformly through the interstices, but forms bubbles which 
push violently through the animated particles so that the bed has the 
appearance of a boiling liquid. With a much greater airflow, whole blocks 
of particles are carried bodily up the tube to several times the original bed 
height, and then disintegrate and fall back, giving a pulsating motion. 
Finally, the mean air velocity in the tube exceeds the velocity of free fall 
of the particles, which are then carried out of the tube. 

A fluidized bed is a very attractive means of bringing a gas into contact 
with a solid, and has lateiy been applied to many chemical processes in which 
contact between a gas and a fixed bed of particles was formerly used. If 
heat is produced by chemical reaction within the bed, temperature gradients 
are reduced by the boiling action of the particles. ‘his means that heat 
can be removed by cooling coils within the bed, or by pneumatically 
conveying the particles to external coolers. When the particles are a 
catalyst, and are subject to fouling, pneumatic conveying can be used to 
remove them from the reacting bed for regeneration. Because of these 
advantages, fluidized beds have tome into use on a very large scale for 
petroleum processing, roasting of ores, and drying of solids. 

This book is a series of chapters on fluidization by leading American 
chemical engineers. ‘The first three chapters deal with fluid mechanics 
and heat transfer characteristics, and the remaining six describe various 
processes using fluidization and also outline some operating and control 
techniques. ‘The first part of the book is therefore of direct interest to 
workers in fluid mechanics. Moreover, they have excellent reasons for not 
ignoring the applications described in the second part; it is good for 
them to know what questions the engineer would really like to answer, 
since even mathematicians usually like to solve useful problems. 
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The first two chapters describe measurements that have been made on 
small-scale fluidized beds and their correlation. Of these correlations, 
the simplest are those which give overall pressure drop and bed height as 
functions of airflow, and have the form shown in figure 1. 4 is the first 
point at which the pressure drop is sufficient to support the weight of the 
bed. ‘Thereafter the bed expands with increasing airflow, but is quiescent 
up to point B. Higher flows induce bubbling and oscillatory motion, and the 
bed height fluctuates between limits shown by the dotted lines. ‘here arises 
the question, are the particles fixed in position between 4 and B? If so, it can 
be presumed that the bed expands by changing from a closely packed con- 
figuration toaloosely packed one. ‘The slight pressure rise from Ato C might 
be caused by particles bridging the tube in modes of closer packing. Point 
B would correspond to the loosest possible packing for particles in contact. 
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Figure 1. Figure 2. 


The characteristics of a particle bed fluidized by water are somewhat 
different (figure 2). In this case there is a unique bed height for all flow 
rates, there are no bubbles, and the particles are evenly spaced. Under 
what circumstances do increments of flow lead to the formation of bubbles 
rather than increased spacing of all the particles? This might be considered 
as a stability problem. Suppose that a number of equal spheres are held 
at points on a uniform lattice, and that fluid is blown through the interstices 
at a rate sufficient to support the weight of each sphere. If a particular 
sphere is released, so as to be free under the action of drag and gravity 
forces, under what conditions will it be stable in the original position ? 
There is experimental evidence to suggest that the relative densities of the 
spheres and fluid would be important, because lead spheres give bubbling 
when fluidized by water, but glass spheres do not; glass beads give bubbling 
when fluidized by air, paper cubes do not! 

The authors of the present book do not give much attention to this 
kind of problem. As engineers, they are concerned with correlations 
which can be used for design, with mechanical features of plant, operational 
details, and applications to various processes. For those who are interested 
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in finding out how fluidization works, the first three chapters are the most 
stimulating. ‘These chapters deal respectively with fluid mechanics, heat 
and mass transfer, and graphical methods of analysis for fluidized systems. 
The first two chapters describe correlations devised in an attempt to relate 
the dimensionless parameters of the fluidized system. ‘This method of 
correlation has been very successful for relatively simple mechanical systems ; 
but for complex systems such as fluid flow through a fixed porous medium, 
its application has been less successful. For example, when matter diffuses 
from the surface of the particles to the fluidizing gas, a mass transfer factor 
j, can be correlated with a modified fluid Reynolds number R. The 
correlation given in chapter 2 covers values of R between | and 10?, and of 
J, between 0-01 and 1-5. Results due to many experimenters are said to 
fall on a single line, although, at a given R, the largest 7, may be six times 
the least. Nevertheless the correlation is useful, and it is difficult to see 
how else to deal with such a complex problem. 

‘The third is perhaps the most original chapter, and deals with a wider 
range of problems than is normally implied by the word fluidization. All 
types of combined fluid and particle flow in a vertical pipe are represented 
on a diagram, with gas velocity as abscissa and vertical pressure gradient as 
ordinate. Roughly speaking, the particles will behave in much the same 
way as a liquid. ‘hus, if they are restrained at the bottom, the particles 
will pack together, and, in the presence of an upward flow of gas, can be trans- 
ferred downwards through a valve, or will How from one limb of a U-tube to 
another under the influence of gravity. As described above, the gas can 
bubble through a fluidized bed of particles, or, if the containing tube is small, 
alternate slugs of gas and particles will be formed. Under different 
conditions, the fluidized bed of particles is *‘ atomized ’, and the separate 
particles travel through the tube at some distance from one another. ‘The 
particles will move upwards in this manner if the gas velocity is very high, 
and downwards if the gas velocity is very low and the particles are not 
restrained at the bottom of the tube. Such an analogy with the behaviour 
of liquids is difficult to express quantitatively, but is useful in designing 
apparatus tor particle transfer. 

The remaining chapters deal with process operation, but are not without 
interest for research workers who wish to come to grips with the real practical 
problems. ‘These chapters would have been better if more numerical 
information about large-scale plant had been given. ‘lhe application to 
large scale plant of the correlations in the earlier chapters is somewhat 
dubious without direct verification from big equipment. 

A general criticism is that the book suffers as the result of being based 


rather too directly on a symposium. ‘Thus each chapter is written in the 
conversational style of its author; the writing is often ungrammatical, and 
many passages are scarcely intelligible, although the colloquial style is 
sometimes refreshing. ‘The book will not become a classic on the subject ; 
it is a progress report which should stimulate the chemical engineer to make 
wider applications, and the research worker to make further efforts to 


understand the mechanism of fluidization. 
J. F. Davipson 

















